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Abstract 

Peridynamics  is  a  nonlocal  formulation  of  continuum  mechanics  in  which  forces 
are  calculated  as  integral  functions  of  displacement  fields  rather  than  spatial  deriva¬ 
tives.  The  peridynamic  model  has  major  advantages  over  classical  continuum 
mechanics  when  displacements  are  discontinuous,  such  as  in  the  case  of  mate¬ 
rial  failure.  While  multiple  peridynamic  material  models  capture  the  behavior  of 
solid  materials,  not  all  structures  are  conveniently  analyzed  as  solids.  Finite  Ele¬ 
ment  Analysis  often  uses  ID  and  2D  elements  to  model  thin  features  that  would 
otherwise  require  a  great  number  of  3D  elements,  but  peridynamic  thin  features 
remain  underdeveloped  despite  great  interest  in  the  engineering  community.  This 
work  develops  nonordinary  state-based  peridynamic  models  for  the  simulation  of 
thin  features.  Beginning  from  an  example  nonordinary  state-based  model,  lower 
dimensional  peridynamic  models  of  plates,  beams,  and  shells  are  developed  and 
validated  against  classical  models.  These  peridynamic  models  are  extended  to 
incorporate  brittle  and  plastic  material  failure,  compounding  the  peridynamic  ad¬ 
vantages  of  discontinuity  handling  with  the  computational  simplicity  of  reduced- 
dimension  features.  These  models  will  allow  peridynamic  modeling  of  complex 
structures  such  as  aircraft  skin  that  may  experience  damage  from  internal  forces 
or  external  impacts. 
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Chapter  1:  INTRODUCTION 


The  final  goal  of  many  mechanical  engineering  analyses  is  the  prediction  and  description  of  ma¬ 
terial  failure.  Even  when  some  amount  of  material  failure  is  acceptable,  conservative  engineering 
models  can  indicate  safe  conditions  by  which  failure  may  be  avoided.  Within  these  operation 
envelopes,  a  wide  variety  of  well-developed  analysis  techniques  are  available  to  predict  mate¬ 
rial  behavior.  Outside  the  envelope,  the  onset  of  failure  may  be  predicted  and  replacement/repair 
indicated,  but  the  actual  progression  of  material  failure  is  unknown.  Because  these  envelopes  con¬ 
servatively  restrict  operation,  absolute  avoidance  of  material  failure  may  require  tradeoffs  such 
as:  reduced  operational  life,  expensive  inspection  and  repair  regimes,  increased  down-time,  and 
reduced  performance  in  other  areas.  Material  models  that  accurately  predict  failure  progression 
can  extend  the  operational  envelope  without  reducing  reliability.  Other  problems,  such  as  those 
related  to  impact,  penetration,  and  blast  resistance,  necessarily  involve  material  failure  and  cannot 
be  accurately  modeled  without  simulating  failure  progression.  Without  reliable  means  for  model¬ 
ing  failure  progression,  these  problems  can  be  tackled  only  by  means  of  extensive  (and  expensive) 
testing  programs.  While  some  extrapolation  from  test  results  is  possible,  many  conditions  of  inter¬ 
est  are  difficult,  expensive,  and/or  dangerous  to  create  and  to  observe.  For  these  reasons,  accurate 
and  reliable  models  for  the  progression  of  failure  for  various  materials  and  conditions  has  long 
been  and  remains  a  major  focus  of  engineering  research. 

1.1  Scope 

This  dissertation  presents  a  nonlocal  equivalent  to  an  Euler-Bemoulli  beam,  along  with  a  method¬ 
ology  for  representing  non-uniform  cross-sections,  plastic  behavior,  and  failure.  Unlike  many  con¬ 
tinuum  beam  theories  that  derive  new  equations  of  motion  (such  as  fourth  order  PDE’s)  from  the 
three-dimensional  elastic  constitutive  model,  or  even  introduce  new  equations  of  motion  without 
reference  to  any  solid  constitutive  model,  the  new  model  is  not  derived  from  prior  ordinary  peridy- 
namic  models  based  on  bond  extension,  but  is  a  material  model  that  directly  resists  bending  defor- 
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mation  while  maintaining  the  same  conservation  of  momentum  equation  as  the  three-dimensional 
model.  In  other  words,  bending  is  a  fundamental  mode  of  deformation  to  this  model,  which  resists 
it  at  the  constitutive  level.  This  beam  model  is  demonstrated  to  be  equivalent  to  Eringen’s  nonlocal 
elasticity  for  small  peridynamic  horizons. 

The  one-dimensional  beam  model  is  then  extended  to  two  dimensions  to  model  the  bending 
behavior  of  flat  Kirchhoff-Love  plates.  The  resulting  1 -parameter  model  is  constrained  to  a  Poisson 
ratio  V  =  1/3.  By  introducing  an  isotropic  bending-state,  the  model  is  extended  to  any  valid 
Poisson  ratio.  The  model  is  combined  with  an  extension-based  model  to  capture  the  effect  of  in¬ 
plane  forces  on  bending  behavior,  resulting  in  a  flat  shell  model  that  is  easily  discretized  to  model 
simple  shapes.  Introducing  the  concept  of  virtual  points  results  in  a  discretized  model  that  is  more 
practically  applicable  and  able  to  model  non-uniform  distributions  of  peridynamic  points,  as  might 
result  from  a  meshing  software.  Using  virtual  points  also  allows  the  simulation  of  curved  shells, 
greatly  expanding  the  class  of  problems  that  can  be  approached  with  this  model.  Because  many 
analyses  of  interest  are  partly  or  wholly  comprised  of  these  types  of  features,  their  development  is 
an  important  addition  to  the  capabilities  of  peridynamic  analysis. 

1.2  Outline 

Chapter  2  of  this  dissertation  reviews  the  literature  for  material  modeling,  focussing  on  the  most 
prominent  PDE-based  material  failure  modeling  techniques,  alternative  nonlocal  models,  peri- 
dynamics,  and  thin  feature  modeling.  Chapter  3  gives  a  short  background  on  Euler  beams  and 
Kirchhoff  plates.  Chapter  4  provides  necessary  background  information  on  peridynamic  models. 
Chapter  5  develops  the  new  bond-pair  and  bond-multiple  peridynamic  bending  models.  Chap¬ 
ter  6  translates  the  continuum  model  into  a  useful  discrete  form  and  extends  that  discrete  model  to 
handle  curved  and  irregular  shapes.  The  results  of  these  discrete  models  are  compared  to  analyt¬ 
ical  and  classical  solutions  for  simple  models  in  section  6.6.  Conclusions  and  avenues  for  future 
development  are  laid  out  in  chapter  7. 
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Chapter  2:  LITERATURE  REVIEW 


2.1  PDE-Based  Failure  Modeling 

Because  material  failure  is  an  important  part  of  many  problems,  there  are  several  computational 
methods  that  attempt  to  simulate  progressing  material  failure.  Most  of  these  methods  start  with 
the  partial  differential  equations  for  conservation  of  momentum  found  in  classical  continuum  me¬ 
chanics,  then  extend  those  models  to  account  for  failure  situations  that  are  poorly  modeled  by  the 
classical  equations. 

One  of  the  most  popular  is  the  extended  Finite  Element  Method  (XFEM),  an  excellent  overview 
of  which  can  be  found  in  [28]  by  Fries  and  Belytschko.  Traditional  finite  elements  use  polynomial 
basis  functions  to  describe  unknown  fields  over  an  element.  As  a  result,  fields  such  as  displace¬ 
ment  always  vary  continuously  across  an  element;  discontinuities  within  an  element  are  poorly 
approximated  even  by  high-order  polynomials.  When  they  are  allowed  at  all,  discontinuities  can 
only  arise  between  elements,  and  are  therefore  very  mesh-dependent;  worse,  to  model  a  moving  or 
evolving  discontinuity  requires  continual  remeshing  around  the  advancing  crack  tip,  adding  greatly 
to  the  computational  complexity  of  the  problem.  The  finite  elements  in  an  XFEM  model  are  “ex¬ 
tended”  or  enriched  by  adding  basis  functions  that  are  discontinuous  or  that  are  continuous  but 
have  discontinuous  derivatives.  In  the  case  of  crack  growth,  this  extension  comes  in  the  form  of 
discontinuous  basis  functions  that  reflect  the  discontinuous  displacement  field  at  a  crack.  In  1982, 
Arnold  [2]  developed  a  mathematical  basis  for  using  discontinuous  finite  elements  and  an  interior 
penalty  method.  This  method  allows  modeling  problems  with  very  high  gradients  (as  at  bound¬ 
aries)  and  greatly  reduces  the  importance  of  meshing  and  the  need  for  remeshing.  While  initially 
demonstrated  for  high-gradient  heat  and  fluid  flow  problems,  the  method  of  adding  discontinuities 
within  an  element  had  strong  potential  for  failure  problems  resulting  in  discontinuous  displace¬ 
ments.  Instead  of  remeshing  along  an  expected  crack  path,  a  process  that  would  require  continual 
remeshing,  the  predicted  direction  of  crack  growth  is  used  to  choose  discontinuous  basis  functions 
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consistent  with  a  crack  along  that  path. 

The  use  of  discontinuous  basis  functions  (such  as  fig.  2.1)  to  model  crack  progression  through 
an  element  was  first  developed  by  Belytschko  and  Black  in  [6]  to  model  straight  and  gently  curving 
cracks  with  less  remeshing.  The  discontinuous  function  in  fig.  2.1  will  look  familiar.  As  the 
transverse  displacement  field  around  the  tip  of  a  Mode  I  crack, 


v{x,y) 


Ki  ry~ 

7^  A  /  cos 
2ii  Y  2// 


K  -  1  +  2  sin^ 


(2.1) 


eq.  (2.1)  and  the  other  crack-tip  displacement  fields  are  a  natural  choice  of  basis  functions.  Longer 
cracks  would  require  remeshing,  but  at  the  root  of  the  crack  rather  than  the  tip.  Root  remeshing 
represents  a  major  improvement  over  tip  remeshing  because  it  needs  to  be  performed  far  less  fre¬ 
quently.  Additionally,  the  path  of  the  advancing  crack  can  travel  at  any  angle  rather  than  following 
the  angle  of  a  mesh  edge.  The  practice  was  then  extended  by  Moes  et  al.  in  [49]  to  use  different 


Figure  2.1:  XFEM  includes  discontinuous  enrichment  basis  functions 

discontinuous  basis  functions  for  cracks  and  crack  tips.  Adding  a  Haar  function  to  model  dis¬ 
continuous  displacements  far  from  the  crack  tip  allowed  cracks  to  pass  through  several  elements, 
enabling  the  modeling  of  cracks  with  greater  length  and  curvature  without  remeshing. 
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As  with  traditional  finite  element  methods,  XFEM  is  mesh-based  and  uses  traditional  finite 
element  basis  functions.  The  elements  are  enriched  extrinsically  by  adding  discontinuous  basis 
functions  to  the  existing  basis  functions  using  the  Partition  of  Unity  concept  described  by  Melenk 
and  Babtiska  in  [44].  The  enrichment  functions  can  capture  arbitrary  discontinuities  in  both  the  pa¬ 
rameter  being  modeled  and  in  its  gradient,  but  the  location  of  the  discontinuities  must  be  predicted 
so  that  the  appropriate  enrichment  functions  can  be  used.  In  [49],  crack  direction  was  predicted 
by  locating  the  maximum  circumferential  stress.  There  are  other  options  for  predicting  a  crack’s 
direction  of  travel,  but  they  do  not  generally  predict  the  formation  of  a  crack  where  none  existed 
before.  Elements  are  enriched  locally  as  illustrated  in  fig.  2.2  rather  than  globally,  to  capture  local 
phenomena  like  crack  growth.  In  some  elements,  every  node  is  enriched,  while  in  others,  only 
some  nodes  are  enriched.  Some  choices  of  enrichment  functions  cause  difficulties  in  the  partially 
enriched  “blending”  elements,  a  problem  for  which  a  variety  of  fixes  are  available  depending  on 
the  situation. 


Figure  2.2:  XFEM  models  enrich  nodes  around  a  discontinuity  [28] 

More  recently,  Sauerland  and  Fries  applied  XFEM  to  two  phase  flow  problems  [62],  including 
standard  dam-breaking  and  rising  droplet  problems.  Holl  et  al.  [32]  used  XFEM  in  the  multi 
scale  modeling  of  crack  propagation,  including  multiple  interacting  cracks.  Mohammadnejad  and 
Khoei  [50]  and  Hunsweck  et  al.  [33]  modeled  using  XFEM  the  combination  of  fluid  and  fracture 
behavior  found  in  hydraulic  fracturing,  a  topic  of  considerable  recent  interest. 

A  second  common  method  of  modeling  material  failure  is  the  Reproducing  Kernel  Particle 
Method  (RKPM),  developed  by  Liu  et  al.  [40].  Mesh-based  models  track  the  connections  between 
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points  in  a  deforming  material,  which  can  become  problematic  when  large  deformations  deform  the 
mesh  and  even  cause  it  to  intersect  itself.  RKPM  is  a  “mesh-free”  method  that  tracks  material  prop¬ 
erties  by  their  values  at  selected  points.  Mesh-free  methods  are  advantageous  for  modeling  large 
deformation  in  fluids  and  solids  because  they  do  not  track  connections  between  points.  Between 
material  points,  properties  are  interpolated  from  their  values  at  nearby  points  by  way  of  integra¬ 
tion  with  a  kernel  function.  The  difference  in  domains  between  RKPM  and  XFEM  can  be  seen  in 
fig.  2.3.  The  key  to  the  RKPM  lies  in  this  kernel  function;  by  using  window  functions  from  wavelet 


Figure  2.3:  Comparison  of  domains  of  influence  for  a)FEM  and  b)RKPM,  by  Guan  et  al.  [31] 


analysis,  a  “reproducing”  kernel  guarantees  that  integrals  of  interpolated  properties  reproduce  the 
integral  of  the  continuous  property  field.  The  “reproducing”  kernel  makes  the  RKPM  a  Partition  of 
Unity  method  like  XFEM,  and  is  a  major  advantage  of  RKPM  relative  to  other  particle  methods. 
It  is  computationally  expensive  to  perform  the  necessary  shape-function  integration,  however.  The 
integration  is  typically  performed  on  a  background  grid,  raising  concern  over  whether  it  is  truly 
“mesh-free”. 

In  its  original  formulation,  the  RKPM  successfully  handles  large  deformations  that  would  cause 
unacceptable  distortion  in  Finite  Element  models.  A  semi-Lagrangian  version  implemented  by 
Guan  et  al.  in  [30]  allows  for  recalculation  of  the  support  function.  This  allows  damage  in  the  form 
of  severing  the  relationship  between  points  that  have  been  pulled  too  far  apart.  Recent  applications 
of  the  RKPM  include  the  work  of  Guan  et  al.  on  fragment-impact  problems  [31]  and  analysis  of 
non-linear  wave  equations  by  Cheng  and  Kim  Meow  [13].  The  RKPM  has  also  been  used  by  Xie 
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and  Wang  [82]  to  analyze  coupled  hydro-mechanical  behavior.  Wu  et  al.  have  coupled  RKPM  and 
FEM  in  their  recent  work  on  fragmentation  and  debris  evolution  [81]. 

Other  methods  of  modeling  failure  include  multiple  finite  element  techniques  and  particle 
methods.  Finite  element  methods  incorporating  “element  death”  remove  from  consideration  el¬ 
ements  that  meet  particular  criteria  and  are  perhaps  the  simplest  way  of  modeling  material  failure. 

Cohesive  zone  elements,  proposed  by  Ortiz  et  al.  [56],  use  element  level  information  to  detect 
the  onset  of  plasticity  (material  instability)  and  add  suitable  deformation  modes  to  model  shear 
banding.  The  additional  deformation  modes  allow  cohesive  zone  elements  to  capture  the  more 
complex  displacements  associated  with  plastic  shear  bands.  Other  cohesive  zone  elements  devel¬ 
oped  by  Needleman  [53]  model  crack  growth,  but  Foulk  et  al.  [27]  note  that  this  requires  very 
fine  meshes  or  prior  knowledge  of  the  crack  path.  They  are  useful  for  cases  such  as  composite 
delamination  or  debonding,  where  the  crack  follows  a  known  surface.  Fang  et  al.  propose  in  [23] 
a  method  of  augmenting  the  cohesive  zone  model  to  work  in  concert  with  XFEM-type  elements 
to  model  both  arbitrary  and  known  crack  paths.  More  recently,  McGarry  et  al.  [43]  and  Mairtm 
et  al.  [41]  developed  a  cohesive  zone  model  to  account  for  crack  closure,  including  crack  surface 
tractions. 

Particle  methods  of  modeling  failure  include  the  Smooth  Particle  Hydrodynamic  (SPH)  method, 
reviewed  by  Monaghan  in  [51],  in  which  a  kernel  (commonly  a  cubic  spline)  is  used  to  create  a 
smooth  interpolation  of  actual  quantities.  Unlike  the  wavelet  basis  functions  of  RKPM,  which 
can  be  made  to  reproduce  a  polynomial  field  of  any  order,  cubic  splines  only  perfectly  reproduce 
constant  fields.  Developed  for  fluids,  Springel  notes  in  [71]  that  it  is  often  used  in  astrophysics 
problems,  where  many  fluid  problems  are  encountered  and  even  “solid"  bodies  deform  under  their 
own  gravity.  It  can  also  predict  elastic  behavior  and  has  been  extended  with  failure  models  by 
Benz  and  Asphaug  in  [7]  by  adding  an  evolving  damage  parameter.  Particle  methods  like  SPH 
handle  fragmentation  very  well,  and  are  used  in  a  variety  of  problems  with  material  interfaces, 
high  strains,  and  multiphase  and  multi-physics  aspects.  Because  cubic  spline  kernels  do  not  re¬ 
produce  a  Partition  of  Unity,  property  interpolation  does  not  accurately  reproduce  the  continuum 
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field.  Additionally,  elastic  SPH  models  do  not  conserve  angular  momentum,  and  in  cases  require 
artificial  viscosities  or  artificial  stresses  to  avoid  numerical  instability. 

The  Material  Point  Method,  developed  by  Sulsky  and  Schreyer  in  [72],  tracks  points  in  both 
Eulerian  and  Lagrangian  meshes.  The  advantage  of  using  both  meshes  is  an  ability  to  handle 
obstacles  and  boundary  conditions  without  difficulty  at  large  deformations.  Wieckowski  notes 
that  the  downside  to  using  both  meshes  is  continual  remeshing  and  added  computational  expense 
[80].  Despite  this  expense,  the  material  point  method  was  extended  to  model  cracks  by  Nairn 
in  [52],  then  further  enhanced  by  Sadeghirad  et  al.  in  [60]  to  improve  stability  when  modeling 
massive  deformations.  The  Lattice  Discrete  Particle  Method  was  developed  by  Cusatis  et  al.  [14, 
15]  to  model  concrete.  Unlike  perviously  mentioned  particle  methods,  the  particles  in  LDPM 
represent  actual  particles  of  aggregate  and  the  cement  between  them,  with  volume  as  well  as  mass. 
It  fills  a  volume  with  variously-sized  particles  generated  from  a  probability  density  function  based 
on  aggregate  size.  The  relationships  between  these  particles  form  tetrahedrons  (the  lattice)  that 
fill  the  volume  and  allow  multi-particle  interaction.  As  with  finite  element  models,  displacement 
between  particle  centers  is  linearized  to  compute  strains,  stresses,  and  forces.  Failure  occurs  at 
predefined  surfaces  between  particles  and  can  include  various  failure  modes.  The  LDPM  is  capable 
of  accurately  modeling  many  failure  conditions  in  concrete,  including  the  impact  of  specimen  size 
on  effective  material  behavior.  These  are  only  a  few  of  the  myriad  of  models  that  attempt  to  model 
the  progression  of  material  failure. 

All  of  these  methods  are  used  to  solve  a  partial  differential  equation  for  conservation  of  mo¬ 
mentum  in  a  continuum.  Because  they  are  based  on  continuum  PDEs,  they  do  not  naturally  develop 
discontinuous  displacements  such  as  cracks.  The  PDEs  that  govern  these  methods  are  ill-defined  at 
the  surface  of  a  crack,  so  cracks  must  be  inserted  within  or  between  elements  after  discretization. 
This  results  in  crack  propagation  that  is  discretization-dependent  as  well  as  computationally  ex¬ 
pensive  and  potentially  unstable.  Although  progress  in  addressing  these  issues  continues,  much  of 
the  difficulty  is  essentially  tied  to  the  undefined  nature  of  derivative  equations  at  discontinuous  dis¬ 
placements.  By  abandoning  the  use  of  displacement  derivatives,  peridynamics  offers  an  alternative 
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way  to  address  discontinuous  displacements. 


2.2  Peridynamic  Modeling 

The  term  peridynamic  was  coined  by  Silling  to  describe  the  new  formulation  of  continuum  me¬ 
chanics  he  developed  in  [64],  From  the  Greek  roots  peri  and  dyna  meaning  near  and  force  re¬ 
spectively,  it  alludes  to  the  nonlocal  force  exerted  by  nearby  points.  In  contrast  to  the  nonlocal 
continuum  mechanics  models  of  Kroner,  Eringen  and  Edelen  [21, 22, 36],  which  formulate  behav¬ 
ior  as  an  integral  function  of  strain  (itself  a  spatial  derivative  of  displacement),  the  peridynamic 
model  casts  material  behavior  at  a  point  as  an  integral  equation  of  the  surrounding  displacement 
rather  than  the  classical  differential  equation.  Because  peridynamic  models  do  not  include  spa¬ 
tial  displacement  derivatives,  discontinuous  displacements  can  arise  naturally  and  can  be  analyzed 
without  first  discretizing  the  problem  or  applying  special  heuristics.  In  classical  continuum  me¬ 
chanics,  linear  momentum  is  conserved  according  to  the  local  eq.  (2.2), 

p(x)u(x)  =  V  ■  P(x)  b(x) ,  (2.2) 

in  which  p  is  the  density,  ii  is  the  second  time  derivative  of  displacement ,  P  is  the  transpose  of 
the  First  Piola  Kirchhoff  stress  tensor,  and  b  is  the  body  force  density,  all  of  which  are  functions 
of  position  x  and  of  time.  Because  P  is  defined  in  terms  of  the  deformation  gradient,  it  is  clear 
that  eq.  (2.2)  is  undefined  for  discontinuous  displacements.  In  fact,  traditional  models  require  even 
the  first  spatial  derivative  of  displacement  to  be  continuous.  Strongly  nonlocal  models  (including 
peridynamics)  replace  the  divergence-of-stress  term  with  an  integral  functional, 

p(x)u(x)  =  f  f(x,q)dKj  +  b(x),  (2.3) 

Jn 

so  that,  instead  of  the  divergence  of  stress,  we  have  the  integral  of  a  “force"  function  f  of  the 
position  vector  x  and  the  position  vector  q  of  a  point  within  the  body  domain  Q.  This  force 
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function  may  depend  on  x,  q,  their  deformed  positions,  the  original  and  deformed  positions  of 
other  points  in  il,  history,  etc.  It  is  common  for  f  to  be  defined  as  0  for  any  pair  of  points  initially 
further  than  S  apart.  The  points  within  5  of  a  point  x  are  the  neighborhood  of  x  and  are  denoted 
in  fig.  2.4  by  %.  By  including  the  behavior  of  nearby  material,  these  models  introduce  an  inherent 


length  scale  to  the  model.  This  length  scale  is  theoretically  determined  by  material  properties, 
though  choice  of  length  scale  is  sometimes  limited  by  computational  demands. 

Constitutive  modeling  of  a  wide  variety  of  materials  is  accomplished  by  choosing  the  appro¬ 
priate  form  for  the  force  function.  The  form  of  the  simplest  such  function  is  a  peridynamic  “bond” 
between  two  points  that  is  modeled  by  a  pairwise  force  function.  While  the  simplest  force  func¬ 
tions  recreate  a  one-parameter  linear  elastic  solid  material,  other  force  functions  can  be  used  to 
model  a  wide  variety  of  material  behaviors,  some  of  which  will  be  outlined  here.  Most  simulation 
of  material  behavior  uses  an  equation  of  motion  reformulated  for  a  discretized  model.  The  typical 
discretization  is  a  mesh-free  numerical  method  in  which  there  are  no  geometrical  connectivities 
between  various  nodes. 

A  force  function  can  be  restricted  to  being  pairwise  (depending  solely  on  the  displacement 
of  the  two  points  x  and  q),  and  still  model  complex  and  varied  behavior.  By  including  a  dam¬ 
age  parameter  that  sets  the  force  contribution  of  “damaged”  bonds  to  0,  Silling  and  Askari  [68] 
were  able  to  model  a  brittle  material  with  natural  crack  formation,  propagation,  and  branching. 
Other  examples  of  damage  propagation  include  impacts  against  brittle  structures  as  in  fig.  2.5 
modeled  by  Demmie  and  Silling  [17]  and  fracturing  of  thermally-stressed  glass  modeled  by  Kilic 

10 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


and  Madenci  [35].  Modeling  progressive  fracture,  including  crack  branching,  is  a  major  advantage 


Figure  2.5:  Peridynamic  model  of  an  airplane  impacting  a  concrete  structure  [17] 

of  peridynamic  formulations.  Using  a  piecewise  force  function,  Dayal  and  Bhattacharya  [16]  were 
able  to  model  phase  transformation  in  ID  and  2D  without  an  additional  constitutive  law;  the  trans¬ 
formations  arose  and  propagated  naturally  as  a  dynamic  instability,  a  result  of  the  force  function 
used.  Peridynamic  models  have  also  been  used  to  analyze  composite  laminates.  In  [83],  Xu  et  al. 
designate  peridynamic  bonds  as  fiber  or  matrix  bonds  with  different  force  functions  to  model  dam¬ 
age  in  composite  laminates.  Kilic  et  al.  model  fiber,  matrix,  and  interfacial  bonds  in  [34]  to  capture 
stacking  order  effects  on  damage  propagation.  Bobaru  [8]  applied  the  peridynamic  model  to  nano 
fiber  networks,  at  a  scale  where  long-range  forces  are  very  apparent.  In  the  same  paper  he  created 
a  Representative  Volume  Element  (RVE)  for  random  networks  of  nano  fibers,  laying  the  ground 
work  for  peridynamic  multi-scale  modeling.  Also  related  to  multi-scale  peridynamic  modeling 
is  work  by  Silling  on  model  coarsening  [65],  fig.  2.6.  An  example  of  a  multi-scale  peridynamic 
simulation  can  be  found  in  [3],  by  Askari  et  al. 

Concrete  is  a  nearly  standard  example  material  in  which  nonlocal  behavior  is  easily  observed, 
and  modeling  the  damage  accumulation  and  proceeding  discontinuity  propagation  has  long  been 
the  goal  of  nonlocal  models  developed  by  Bazant  and  Pijaudier-Cabot  [5]  among  others,  signifi¬ 
cantly  predating  peridynamics.  In  [29],  Gerstle  et  al.  use  rotational  degrees  of  freedom  to  create 
a  concrete  material  model,  capable  of  describing  a  linear  elastic  material  with  any  Poisson  ratio, 
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Figure  2.6:  Silling’s  illustration  of  course-graining  in  time  from  [65]. 

that  also  handles  material  failure.  Peridynamic  models  are  not  limited  to  force-displacement  re¬ 
lationships;  the  theory  has  also  been  applied  to  diffusion  processes  and  multiphysics  problems. 
Peridynamic  models  can  simulate  heat  transfer  [9]  and  diffusion  [10]. 

Mathematical  analyses  of  simplified  cases  have  also  been  fruitful.  Weckner  [78]  determined 
analytical  solutions  to  the  infinite  bar  problem.  Emmerlick  and  Zimmerman  proved  solution  ex¬ 
istence  and  uniqueness  in  the  simplest  case  of  the  peridynamic  bar  [20].  Mikata  found  additional 
analytical  solutions  for  the  bar  problem  [46].  In  3D,  Weckner  constructed  Green’s  functions  for  an 
infinite  peridynamic  solid  in  [79].  All  of  this  work  was  done  with  peridynamic  models  limited  to 
pairwise  force  functions. 

Other  than  Gerstle’s  aforementioned  micropolar  peridynamic  model,  the  pairwise  force  func¬ 
tion  limits  3D  solid  materials  to  a  Poisson  ratio  of  u  =  1/4.  To  model  additional  material  behavior. 
Silling  et  al.  generalized  the  underlying  peridynamic  concept  of  bonds  and  forces  and  introduced 
state-based  peridynamic  models  in  [66].  By  freeing  the  force  function  from  the  pairwise  restric¬ 
tion,  state-based  models  allow  the  force  relationship  between  two  points  to  depend  on  the  collec¬ 
tive  behavior  of  all  nearby  material.  Using  the  concept  of  a  deformation  vector-state  allows  for  the 
construction  of  correspondence  models  that  can  recreate  any  classical  constitutive  model.  These 
correspondence  models  use  the  deformation  state  to  approximate  the  deformation  gradient  tensor, 
then  use  the  deformation  gradient  tensor  to  calculate  force  contributions.  State-based  models  were 
used  by  Foster  et  al.  to  simulate  viscoplasticity  and  hardening  in  [26],  and  rate  dependent  failure 
in  [25],  with  others,  via  an  energy  criterion.  Mitchell  describes  state-based  models  for  plasticity 
in  [47]  and  viscoelasticity  in  [48].  A  non-ordinary  state-based  model  was  used  by  Warren  et  al.  to 
simulate  fracture  in  [77].  More  recently,  Tupek  et  al.  have  incorporated  the  idea  of  peridynamic 
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damage  into  a  Johnson-Cook  based  damage  state  that  accumulates  with  plastic  strain  [75]. 


2.3  Other  Nonlocal  Elasticity  Models 

The  peridynamic  formulation  of  continuum  mechanics  is  neither  the  only  nor  the  first  nonlocal 
model.  Nonlocal  elasticity  generally  allows  for  forces  at  a  point  that  are  dependent  on  the  material 
configuration  of  an  entire  body,  rather  than  the  configuration  at  that  point  [22].  While  long-range 
forces  are  obvious  at  the  molecular  model,  material  at  larger  scales  is  conventionally  modeled  as 
though  internal  forces  are  local  or  contact  forces  [36].  The  result  of  such  approximation  is  accu¬ 
rate  for  deformations  that  are  homogeneous,  but  introduces  some  inaccuracy  for  inhomogeneous 
deformations  like  the  propagation  of  waves  with  short  wavelengths.  One  way  to  distinguish  be¬ 
tween  homogeneous  and  inhomogeneous  deformations  is  to  incorporate  higher-order  gradients  of 
deformation.  While  stress  in  classical  elasticity  is  a  function  of  the  (first)  gradient  of  deformation, 
Eringen’s  formulation  of  a  nonlocal  modulus  in  [21]  approximates  a  weighted  sum  of  the  first  and 
second  order  gradients.  This  introduces  a  length  scale  to  the  model  and  has  the  effect  of  smear¬ 
ing  out  local  deformation  inhomogeneities  over  the  surrounding  material,  while  maintaining  the 
conventional  result  for  homogeneous  deformations. 

Previous  work  in  the  nonlocal  mechanics  of  beams  is  motivated  by  the  observed  stiffening  of 
nanoscale  cantilevers.  Challamel  and  Wang  demonstrate  in  [12]  that  Eringen  nonlocal  elasticity 
cannot  reproduce  the  scale  stiffening,  but  that  stiffening  does  result  from  other  gradient-elastic 
models  and  models  incorporating  nonlocal  curvature.  Because  all  of  these  models  incorporate 
higher-order  gradients  of  deformation,  they  impose  stronger  continuity  requirements  than  classical 
elasticity,  and  are  unsuitable  for  discontinuous  displacements.  Because  the  gradients  are  evaluated 
locally,  gradient  models  are  called  weakly  nonlocal.  Recent  work  by  Paolo  et  al.  [57]  develops  a 
displacement-based  beam  in  which  relative  axial  displacement,  shear  displacement,  and  rotation  of 
non-adjacent  beam  segments  are  resisted  by  three  kinds  of  nonlocal  spring,  whose  stiffnesses  can 
be  tuned  to  the  expected  material  behavior.  With  the  appropriate  nonlocal  stiffnesses,  their  model 
reproduces  the  nanoscale  cantilever  stiffening  effect. 

13 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Similarly,  Duan  and  Wang  [19]  applied  Eringen-type  elasticity  to  the  quasi- ID  problem  of 
axisymmetric  bending  in  nanoscale  plates.  Pradhan  and  Murmu  [58]  extended  the  concept  to 
buckling  in  single-layered  graphed  sheets,  a  fully  2D  problem.  Later,  Ansari  et  al.  [1]  modeled  the 
vibration  of  single-layered  graphed  sheets  using  Eringen-type  elasticity. 

Nonlocal  effects  have  also  been  incorporated  into  many  of  the  modeling  techniques  previously 
discussed.  Bazant  and  Chang  incorporated  nonlocal  strain-softening  into  a  finite  element  model 
in  [4].  Any  interpolating  particle  method  will  exhibit  some  measure  of  nonlocality,  but  some  ex¬ 
plicitly  model  nonlocal  phenomena.  Vignjevic  et  al.  used  SPH  to  model  nonlocal  strain- softening 
in  [76],  and  in  [1 1],  Burghardt  et  al.  developed  a  material  point  method  that  incorporates  nonlocal 
plasticity. 

2.4  Thin  Features 

Many  engineering  analyses  concern  shapes  that  have  one  dimension  much  greater  than  another; 
numerical  modeling  the  behavior  of  these  shapes  can  be  a  considerable  challenge  for  methods 
designed  for  3D  solids.  In  finite  element  models,  for  example,  calculations  can  become  unstable  or 
too  stiff  when  individual  elements  become  long  and  thin.  To  avoid  such  elements  while  maintaining 
model  fidelity  requires  a  very  large  number  of  solid  elements.  By  making  some  assumptions 
about  the  behavior  along  the  thin  direction,  many  such  shapes  can  be  modeled  as  ID  beams  or  2D 
plates  or  shells  without  great  loss  of  accuracy.  A  comprehensive  review  of  the  classical  continuum 
mechanics  associated  with  thin  features  by  Reddy  [59]  also  includes  a  section  on  the  finite  element 
analysis  of  plates  and  shells.  Material  failure  in  classical  thin  features  is  modeled  using  the  same 
techniques  as  in  solids.  Dolbow  et  al.  use  XFEM  in  [18]  to  model  fracture  in  plates.  Li  et  al.  use 
a  variant  of  RKPM  in  [37]  to  model  plastic  deformation  in  shells.  More  recently,  Xu  et  al.  have 
applied  XFEM  to  plate  plasticity  problems  [84],  and  Memar  Ardestani  et  al.  have  used  RKPM  to 
model  functionally  graded  plates  [45].  Other  authors  use  cohesive  zone  elements  [38]  or  SPH  [42] 
to  study  failure  in  thin  features. 
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2.4.1  Peridynamic  Models 


Reduced  dimension  thin  features  such  as  bars  [20,46,69,78],  plates  [35],  and  membranes  [70]  have 
been  modeled  using  peridynamics,  but  these  models  are  used  for  in-plane  or  membrane  forces  as 
shown  in  fig.  2.7.  Because  traditional  peridynamic  models  exert  forces  in  the  direction  of  the  bonds 


Figure  2.7:  Tearing  a  peridynamic  membrane  [70] 


between  points,  they  are  not  well-suited  for  bending  problems  of  thin  shapes,  in  which  force  and 
displacement  are  both  nearly  perpendicular  to  bonds  connecting  material  points  at  separate  points 
on  a  surface.  Just  as  with  solid  finite  elements,  most  peridynamic  models  of  thin  features  like  the 
tubes  in  fig.  2.8  have  included  several  nodes  through  the  thickness  of  a  thin  part  to  capture  bending 
behavior.  Also  as  with  solid  finite  elements,  this  leads  to  very  fine  discretization  of  thin  features, 
even  when  the  expected  behavior  is  quite  simple.  This  greatly  increases  the  computational  expense 
of  modeling  parts  with  thin  features. 

A  recent  paper  by  Taylor  and  Steigmann  [73]  partially  addresses  this  issue  by  starting  with 
mathematical  analysis  of  the  continuous  3D  bond-based  peridynamic  solid  model  of  a  thin  plate. 
By  using  a  continuous  model,  they  avoid  the  difficulties  associated  with  discretizing  thin  fea- 
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Before 


After 

(brittle  model) 


(plastic  model) 


Figure  2.8:  A  peridynamic  cylinder  uses  several  nodes  through  its  thickness  in  [39] 


tures.  Applying  asymptotic  analysis  to  the  continuous  model,  they  reduce  a  3D  solid  model  to 
2  dimensions.  The  asymptotic  reduction  is  accomplished  by  the  addition  of  degrees  of  freedom 
for  the  derivative  of  displacement  with  respect  to  the  through  thickness  direction.  By  making 
the  through-thickness  derivative  of  displacement  vector  an  independent  variable,  the  resulting  flat 
model  includes  a  measure  of  angular  deformation  that  allows  resistance  to  bending.  Using  a  sim¬ 
ple  bond-stretch  damage  criterion,  the  Taylor  and  Steigmann’s  reduced  model  was  able  to  capture 
the  out-of-plane  displacement  (fig.  2.9)  associated  with  crack  propagation  behavior  (fig.  2.10)  in 
a  pre-cracked  plate  under  tension  loading.  In  general,  however,  the  asymptotic  reduction  model 
encounters  difficulty  when  nonlinear  behaviors  like  damage  are  implemented.  The  use  of  a  bond- 
stretch  criterion  as  implemented  is  only  appropriate  when  deformation  is  dominated  by  in-plane 
tension,  as  failure  caused  by  bending  will  not  be  captured.  Because  of  its  basis  in  the  3D  bond- 
based  solid  model,  Taylor  and  Steigmann’s  model  is  limited  to  a  Poisson’s  ratio  of  Vs.  While  it 
is  possible  that  future  analysis  will  extend  the  asymptotic  reduction  to  state-based  model,  allow¬ 
ing  for  arbitrary  Poisson’s  ratios,  there  are  significant  mathematical  hurdles  that  will  have  to  be 
overcome. 
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Figure  2.9:  Taylor  and  Steigmann’s  asymptotic  reduction  allows  for  bending  resistance  in  a  2D 
plate  in  tension  [73] 
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Chapter  3:  CLASSICAL  BACKGROUND 


3.1  Euler-Bernoulli  Beam  Theory 


The  simplest  representation  of  bending  behavior  is  found  in  the  Euler-Bernoulli  beam  equation, 


=  Q, 


(3.1) 


which  describes  the  transverse  deflection  of  long-slender  beams.  Equation  (3.1)  is  for  a  beam  in 
the  X  direction,  whose  displacement  v  is  perpendicular  to  x  and  parallel  to  the  distributed  load  q 
acting  on  the  beam.  While  it  is  only  applicable  to  small  deformations  and  rotations,  there  are  many 
problems  that  can  be  easily  and  usefully  simplified  by  the  application  of  Euler  beams. 

It  is  important  to  note  that  eq.  (3.1)  is  not  a  constitutive  equation;  it  does  not  relate  stress  and 
strain  for  a  material.  Instead,  it  is  derived  from  the  constitutive  equation  for  a  linearly-elastic  ma¬ 
terial  and  some  assumptions  about  the  deformation.  Euler-Bernoulli  beam  theory  is  much  simpler 
than  a  3D  analysis  of  a  beam  as  a  solid  material,  but  it  requires  the  following  assumptions  to  be 
met: 


•  Slenderness:  the  length  of  the  beam  should  be  20  times  its  other  dimensions 

•  Loaded  transversely,  no  axial  loads  or  torques 

•  Small  deformations  and  rotations 

•  Plane  cross-sections  of  the  beam  remain  plane  under  deformation 

•  Initially  straight,  and  symmetric  about  the  plane  of  bending 

Eirst,  we  take  a  infinitesimal  slice  of  a  beam  as  depicted  in  fig.  3.1.  Because  we  will  want  to 
compare  to  this  model  later,  we  use  a  different  ^-coordinate  direction  than  many  presentations.  We 
use  this  infinitesimal  slice  to  derive  the  relationship  between  load,  shear,  and  moment. 
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da; 

Figure  3.1:  Infinitesimal  Beam  Slice 

A  force  balance  demonstrates  that  the  x-derivative  of  internal  shear  force  is  the  applied  load. 

V  =  V  +  +  ^  =  -q.  (3.2) 

da; 

The  moment  balance  can  be  performed  about  any  point,  but  it  is  simplest  to  use  the  right  side. 

d-T  diA/[ 

M-Vdx  +  q  dx—  =  M  +  dM  ^  —  =  -(/  +  0(dx) .  (3.3) 

2  da; 

Next,  we  will  use  the  deformation  assumptions  to  frame  the  bending  as  a  simple  arc  from  which 


Figure  3.2:  Small  deformation  in  an  Euler  beam 

we  can  deduce  a  relationship  between  the  transverse  displacement  of  a  beam  section  and  its  angle 
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and  radius  of  curvature  to  be 


1  d(b 


=  K  . 


(3.4) 


p  dx  dx^ 

The  radius  of  curvature,  p  is  the  radius  of  the  arc  formed  by  the  beam’s  neutral  axis,  which  does  not 
change  in  length  as  the  beam  bends.  We  use  this  same  arc  to  find  the  strain  in  material  that  is  at  a 
distance  y  from  the  neutral  axis.  In  a  linearly  elastic  material,  the  resulting  stress  is  proportional  to 
Young’s  modulus  E.  To  find  the  bending  moment  that  results  from  that  stress  profile,  we  integrate 
the  moment  contribution  through  the  beam. 


^top  ^top 

M-esisted  =  /  yc^dA  =  /  End  A  . 

J  bottom  J  bottom 


(3.5) 


In  equilibrium,  the  moment  resisted  will  be  equal  and  opposite  to  the  moment  applied 


K  = 


M 


applied 


E!'Z.V^dA 


(3.6) 


The  integral  in  the  denominator  is  the  bending  resistance  of  the  beam’s  shape,  called  its  second 
moment  of  area,  about  the  neutral  axis.  The  second  moment  of  area  is  generally  represented  by 
/.  By  combining  eqs.  (3.2)  to  (3.4)  and  (3.6),  we  reproduce  eq.  (3.1),  the  Euler-Bernoulli  beam 
equation.  If  the  Young’s  modulus  and  beam  cross  section  are  constant  throughout  the  beam,  the 
equation  can  be  further  simplified  to 


^  dS 

Obviously,  to  solve  for  the  deformed  position  we  will  require  4  constraints  from 
conditions.  The  nature  of  the  constraints  is  determined  by  the  support  configuration, 
common  end  conditions  are  listed  table  3.1.  Other  end  constraints,  such 


boundary 
The  most 
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Table  3.1:  Common  Beam  End  Conditions 


End  Condition 

Contraint  1 

Constraint  2 

End  Load 

II 

v"  =  0 

End  Torque 

u'"  =  0 

II 

-d 

Simply  Supported 

v"  =  0 

^  ^support 

Clamped 

y  =  <lamp 

^  ~  ^clamp 

3.2  Kirchhoff-Love  Plate  Theory 


Like  Euler-Bernoulli  beam  theory,  Kirchhoff-Love  Plate  theory  begins  by  making  simplifying  as¬ 
sumptions  about  the  plate  and  its  deformation. 

•  Slenderness:  the  length  and  width  of  the  plate  should  be  20  times  its  thickness 

•  Loaded  transversely,  no  in-plane  loads  or  torques 

•  Small  deformations  and  rotations 

•  Straight  lines  normal  to  the  plate  remain  normal  to  the  center  plane  of  the  plate  under  defor¬ 
mation 

•  Initially  flat 

These  assumptions  allow  us  to  simplify  our  strain-displacement  relationships  and  formulate  our 
curvatures  and  strains  in  terms  of  the  transverse  displacement  w; 


d^w 


Kini 


1^X11 


^  “  Qy2  ’ 


dxdy  ’ 


Sy  ZKy  ,  ^xy  ‘2zK 


xy  ■ 
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To  find  the  stresses  at  any  point,  we  apply  Hooke’s  law  and  find 


(J7.  —  — 


Gy  — 


Ez 
1  — 
Ez 

1  — 

Ez 


(^Kx  Pt^y')  , 
{Ky  +  PKx)  , 


Txy 


1  +  P 


-K. 


xy 


As  with  the  beam,  these  stresses  can  be  integrated  over  the  thickness  of  the  plate  to  determine  the 
resulting  moments. 


nt/2 

II 

J  -t/2 

My  = 

ntl2 

J  -  t/2 

nt/2 

II 

ZGxdz 


ZGydZ 


ZTxydZ 


E  (^fxx  T  PKy^  , 

-D  {Ky  +  PKx)  , 
E  {i  p)  t^xy  ) 


with 

“  12  (1  -  1/2) 

as  the  flexural  rigidity.  Note  that  a  beam  of  unit  width  would  have  a  second  moment  of  area  I  =  ^, 
so  the  flexural  rigidity  of  a  plate  is  greater  by  a  factor  of  The  difference  between  the  values 
reflects  the  fact  that,  when  a  plate  is  bent  in  one  direction,  poisson  ratio  effects  promote  curvature 
in  the  opposite  direction  at  a  right  angle.  Constraining  this  saddle-shaped  (or  anticlastic)  bending 
increases  the  stiffness  or  rigidity  of  the  plate. 

Constructing  the  governing  equation  is  more  complicated  than  for  beams  because,  instead  of 
shear  V  and  moment  M  on  an  infinitesimal  section,  we  have  shears  I4  and  Vy  and  moments  Mx, 
My,  and  Mxy.  Balancing  forces  gives  us 


dx  dy 
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for  transverse  pressure  p.  Balancing  moments  gives 


about  X  : 
about  y  : 


dM,y  ^  dMy 

dx  dy 

dM^y  ^  dM, 
dy  dx 


Vy  =  0, 
14  =  0. 


Combining  the  force  and  moment  balance  equations  generates  the  governing  equation 


d^w  d^w  p 
dx^  dx‘^dy^  dy^  D 


As  with  the  beam  case,  we  will  need  to  apply  boundary  conditions  in  order  to  find  the  deformed 
configuration. 
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Chapter  4:  PERIDYNAMICS  BACKGROUND 


4.1  Peridynamic  States 


Introduced  by  Silling  et  al.  in  2007  [66],  peridynamic  states  are  functions  of  the  behavior  of  the 
continuum  points  surrounding  each  location.  As  is  appropriate  for  a  theory  based  on  force-carrying 
bonds,  states  often  operate  on  vectors.  The  most  common  states  are  scalar-states  and  vector-states 
which  are  scalar  and  vector  valued,  respectively.  As  a  matter  of  convention,  scalar  states  are  usually 
denoted  by  lowercase  letters  (e.g.  a,  b),  while  vector  states  are  denoted  by  uppercase  letters  (e.g. 
A,  B).  A  state  that  operates  on  vectors  and  is  itself  vector  valued  naturally  brings  to  mind  a 
two-point  tensor  such  as  the  deformation  gradient;  unlike  a  second  order  tensor,  which  can  only 
map  vectors  linearly  to  other  vectors,  vector-states  can  produce  nonlinear,  discontinuous,  or  even 
noninvertable  mappings.  This  difference  is  illustrated  in  fig.  4.1. 


Symmetric  tensor  A 


Figure  4.1:  The  deformation  tensor  linearly  maps  spheres  to  ellipsoids,  while  a  vector  state  can 
map  spheres  nonlinearly  to  complex  and  even  discontinuous  shapes  [66] 


The  mathematical  properties  of  states  and  several  related  operators  are  defined  in  [66].  Im¬ 
portant  properties  of  states  are  magnitude  and  direction,  while  important  operations  include  the 
addition  and  composition  of  states,  inner  and  tensor  products,  and  the  Frechet  derivative  of  a  func- 
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Table  4.1:  Common  State  Operation  Nomenclature 


Operation 

Notation 

Meaning 

Addition 

(a  +  b)  {$) 

+  h.{^) 

Multiplication 

(M)(0 

Scalar  Product 

(A-B)(0 

A(«)  ■  B{«) 

Composition 

(AoB)«) 

A(B(^)) 

Dot  Product 

A.B 

[  A(o  ■  m) 

Jn 

a*h 

Jn 

Vector  Norm 

lAKO 

\m)\ 

State  Norm 

IIAII 

Frechet  Derivative 

V/(A) 

df 

dA 

dA 

tion  with  respect  to  a  state.  While  some  of  these  operations  are  intuitive,  the  nomenclature  may  not 
be.  Refer  to  table  4.1  for  the  notation  of  common  state  operations.  For  a  more  rigorous  definition 
and  examples  of  the  Frechet  derivative,  see  appendix  A. 

4.2  State-based  Models 

Conservation  of  linear  momentum  in  the  state-based  peridynamic  formulation  results  in  the  equa¬ 
tion  of  motion 

p(x)u(x)  =  [  (T[x](q  -  x)  -  T[q](x  -  q))dVq  +  b(x) ,  (4.1) 

JQ, 

in  which  T[  ](  )  is  a  force  vector-state  that  maps  the  vector  in  angle  brackets,  (),  originating  at 
the  point  in  square  brackets,  [  ],  to  a  force  vector  acting  on  that  point.  The  deformed  image  of  the 
vector  (q  —  x)  is  defined  as  the  deformation  vector-state,  usually  denoted  Y  and  formulated 

Y[x]  (q  -  x)  =  (q  -  x)  +  (u(q)  -  u(x))  (4.2) 
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for  a  displacement  field  u.  Just  as  stress  and  strain  are  work  conjugate,  so  too  are  the  force  and 
deformation  vector  states  for  hyperelastic  materials. 

State-based  models  include  surrounding  material  behavior  illustrated  in  fig.  4.2  in  the  force 
function  between  each  pair  of  continuum  points.  It  is  common  for  the  formulation  of  the  force 
state  T  to  be  scaled  by  a  weighting  function,  commonly  represented  by  lo,  that  makes  explicit  the 
region  in  which  the  force  relationship  between  points  is  nonzero.  Perhaps  the  simplest  and  most 
common  weight  function  is  eq.  (5.8),  representing  a  constant  nonzero  value  for  bonds  shorter  than 
the  peridynamic  horizon  5. 

|l 

u($)  =  I  (4.3) 

[o  if\^\>S. 

Many  evaluations  of  forces,  energies,  and  other  states  at  a  given  point  are  reduced  from  being  inte¬ 
grals  over  the  entire  body  to  being  integrals  over  that  point’s  neighborhood  H.  This  is  particularly 
useful  when  trying  to  apply  peridynamic  models  to  real-world  problems,  especially  when  using 
computer  models.  If  the  force  state  T  is  always  in  the  same  direction  as  the  deformation  state 


Figure  4.2:  The  body  (2  deformed  by  the  deformation  state  Y 


Y,  then  the  force  exerted  by  a  “bond”  between  points  is  in  the  same  direction  as  the  deformed 
bond,  and  the  model  is  called  ordinary.  Ordinary  state-based  models  can  reproduce  linear  elastic 
materials  with  arbitrary  Poisson  ratios  by  separating  dilatory  and  deviatoric  deformations  and  the 
energy  corresponding  to  each.  They  can  also  model  a  variety  of  elastic  and  inelastic  behaviors. 

There  is  no  requirement  that  force  states  be  in  the  same  direction  as  their  associated  deforma- 
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tion  states,  and  models  in  which  they  are  not  in  the  same  direction  are  called  nonordinary.  Some 
additional  care  is  needed  to  ensure  that  angular  momentum  is  conserved  in  a  nonordinary  state- 
based  model,  but  they  are  still  perfectly  legitimate.  Silling  et  al.  demonstrates  the  possibility  of 
such  models  in  [67],  but  very  little  work  has  touched  on  their  use.  Foster  et  al.  [26]  and  Warren  et 
al.  [77]  show  that  some  correspondence  models,  which  approximate  the  deformation  gradient  and 
use  it  to  calculate  bond  forces,  result  in  non-ordinary  state-based  constitutive  models  for  finite  de¬ 
formations.  Bond-based,  ordinary  state-based,  and  nonordinary  state-based  models  are  illustrated 
in  fig.  4.3. 


Non-ordinary  state  based 


Ordinary  state  based 

Figure  4.3:  Illustration  of  the  three  types  of  peridynamic  models,  from  specific  to  general  [66] 


It  should  be  clear  that  many  of  the  concepts  of  classical  continuum  mechanics  have  direct 
equivalents  in  peridynamic  modeling.  Table  4.2  lays  out  some  of  the  simplest  parallels  between 
classical  and  peridynamic  formulations. 
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Table  4.2:  Peridynamic  Equivalents  of  Classical  Concepts 


Concept 

Classical 

Peridynamic 

Kinematics 

F 

Y 

Linear  Momentum 

Vo- 

/  (T[x]  (q  -  x)  -  T[q]  (x  -  q))dVq 

JQ 

Angular  Momentum 

T 

a  =  a 

[  Y[x]  (q  -  x)  X  T[x]  (q  -  x)dYq  =  0 
Jn 

Constitutive  Law 

cr  =  <T(e) 

T  =  T(Y) 

Stress  Power 

ecT 

T*  Y 

4.3  Bond-based  peridynamics 

If  the  force  state  T[x](^)  depends  only  on  the  deformed  positions  of  the  points  at  the  end  of  the 
bond  then  the  model  is  called  bond-based.  In  bond-based  peridynamic  models,  each  pair  of 
points  is  treated  separately,  without  consideration  of  the  behavior  of  other  points.  This  makes 
bond-based  models  much  simpler  computationally  than  general  state-based  models,  and  reduces 
the  equation  of  motion  to 


p(x)u(x)  =  [  f(u(q)  -  u(x),  q  -  x)dl/q  -f  b(x) .  (4.4) 

Ja 

By  choosing  an  appropriate  function  f ,  this  model  can  reproduce  the  results  of  linear  elasticity  for 
solid  materials  with  a  Poisson  ration  u  =  1/4  and  2-dimensional  materials  with  a  Poisson  ration 
0  —  1/3.  It  can  also  be  used  to  investigate  a  range  of  nonlinear  behaviors  by  changing  the  force 
function  (examples  in  fig.  4.4).  To  conserve  momentum  in  a  bond-based  model,  it  is  only  necessary 
that  f  satisfy 

f  (u(q)  -  u(x) ,  q  -  x)  =  -f  (u(x)  -  u(q) ,  x  -  q) ,  (4.5) 

i.e.  the  forces  exerted  at  the  opposite  ends  of  the  bond  between  x  and  q  must  be  equal  and  opposite. 
The  first  peridynamic  models  were  all  bond-based,  and  provide  useful  insight  into  many  complex 
material  failure  phenomenon  despite  their  limitations. 
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Figure  4.4:  Bond-based  models  can  describe  a  variety  of  material  behaviors 

4.4  Important  Peridynamic  Models 

Though  they  cover  only  a  small  portion  of  the  behaviors  modeled  with  peridynamics,  these  few 
examples  should  serve  to  illustrate  the  form  and  analysis  of  peridynamic  material  models. 

4.4.1  Bond-based  Elastic  Solid 

The  simplest  peridynamic  model  treats  each  bond  as  a  linear  spring  between  two  points.  In  the 
bond-based  formulation,  there  is  no  interaction  between  different  bonds,  so  the  force  function  is 

f  (u(q)  -  u(x),  q  -  x)  =  w(|q  -  x|)  c  s  [(q  +  u(q))  -  (x  +  u(x))]  ,  (4.6) 

with  weighting  function  uj,  spring  constant  c,  and  the  stretch  s  defined  by 

s  =  |(q-hu(q))  -  (x-Fu(x))|  -  |q  -  x|  .  (4.7) 
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For  a  deformation  gradient  F  representing  small  uniform  displacements,  i.e.  Vu  <<  1,  the  stretch 
of  bond  ^  =  q  —  X  is 

For  reasons  that  will  become  clear  in  the  discussion  of  the  next  model,  we  calibrate  the  spring 
constant  c  following  the  approach  of  [66],  by  comparing  the  energy  to  that  of  a  classical  solid 
under  purely  deviatoric  deformation,  so  that  =  efj.  The  energy  of  this  spring  will  be  in  units  of 
energy  per  volume  squared,  so  that  integration  over  all  the  springs  at  a  point  gives  energy  per  unit 
volume. 


c  s 


w  = 


W  = 


2  ’ 


'n 


^(lei)l^l  dV,, 


—  d 


lei  j  V  lei 


(4.8) 


Because  u  depends  only  on  |^|,  we  can  rewrite  this  integral  in  spherical  coordinates  as 


W  =  -eU 


2'«'“  L 


p27r  P7V 


d  A 


'0  ^0 


sin(0)  dcf)  dO  dr  . 


Recognizing  that  —  r  sin  0  cos  0,^2  =  r  sin  (j)  sin  6*,  ^3  =  r  cos  0,  we  can  see  that  configurations 
of  [i,  j,  k,  1]  with  an  odd  number  of  any  index  result  in  integrals  with  an  odd  number  of  one  or  more 
of  cos  9,  sin  9,  cos  0,  and  therefor  are  equal  to  0.  For  the  remaining  configurations. 


uj{r) 


r*27r  /*7r 


/o  JO 


47r 

(r^  sin^  0  cos^  9  sin^  9)  sin(0)  d9  dr  =  — 

lo 


This  leaves  only  configurations  such  as  [1, 1, 3, 3],  [1, 2, 1, 2]  and  [3, 2, 2, 3], 
by  (SikSji  +  SiiSjk  +  6ijSkl)  .  Additionally,  any  combination  with  i  =  j  or  k 


C(j(r)r^  dr . 


which  we  can  indicate 
=  I  results  in  terms 
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or  efi^.  Such  terms  sum  to  0  in  deviatoric  deformation,  leaving  only  {SikSji  +  SuSjk). 


d  47r 


W  dr{SikSji  +  Sudjk) , 


2  *^'“15  Jo 


dvr 


=  c6j6j-  /  u{rydr, 


15 


ij  ij 

15 


m . 


To  force  the  result  to  be  independent  of  the  horizon  <5,  we  normalize  the  expression  by 


m=  <^(|^|)|^|  =4:77  u>{r)  dr  . 


(4.9) 


Jn  Jo 

By  comparing  to  the  classical  strain  energy  density  Vl  =  ji  for  shear  modulus  /i,  we  can 
determine  the  appropriate  bond  stiffness, 


c  = 


15 /i 
m 


Applying  a  purely  dilational  deformation  to  the  same  model  is  far  easier.  With  dilation  O^d,  the 
stretch  of  a  bond  in  any  direction  is 


s  = 


5 


and  the  corresponding  energy  is 


n2  rS  p27r  ptt 

=  — /  u{r)  /  /  r'^  sm{(p)  d(p  d9  dr  , 

18  Jo  Jo  Jo 


C  (*30 
— —m , 
18 

15 


(4.10) 
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This  shows  that  the  model  based  on  bond-stretch  has  a  bulk  modulus  that  is  of  its  shear  modulus, 


indicating  a  Poisson’s  ratio  of  'A. 

A  nearly  identical  analysis  can  be  performed  on  a  simpler  2D  version  of  the  same  model,  de¬ 
parting  after  eq.  (4.8).  Applying  a  purely  deviatoric  in-plane  shear  to  such  a  model,  and  comparing 
the  resulting  energy  to  that  of  a  classical  plate  with  thickness  t  gives  us 

c=— m2D=  f  . 

"^2d  Jn2D 

Applying  a  planar  dilation  deformation  to  a  2D  plate  results  in  strain  energy  consistent  with  a 
Poisson’s  ratio  of  '/s  rather  than  the  value  of  !4  found  for  the  3D  solid. 

4.4.2  State-based  Elastic  Solid 

The  state-based  linear  isotropic  peridynamic  solid  material  model  is  both  important  and  illustra¬ 
tive.  Developed  in  [66],  it  uses  many  of  the  important  characteristics  of  peridynamic  states  to 
model  a  linearly-elastic  material  with  any  valid  Poisson’s  ratio.  The  extension  state  e  is  exactly 
the  same  as  the  stretch  s  in  eq.  (4.7).  Classical  material  models  dealing  with  metal  plasticity  often 
separate  deformation  into  dilation  and  deviation  components.  Similarly,  the  extension  state  can  be 
decomposed  into  isotropic  and  deviatoric  extension  states.  Using  m  from  eq.  (4.9)  as  above  for 
normalization  the  dilation  state  is  defined 

hD  =  -  i  ^(mkdVi. 

^  Jn 

The  isotropic  and  deviatoric  extension  states  are  defined  in  turn 

3  -  -  - 


(4.11) 


(4.12) 
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If  the  energy  associated  with  dilation  is  set  to 


ly*  = 


k 


then  the  corresponding  force  state  is 


t  = 


SkOsD 

m 


We  saw  in  the  analysis  of  the  bond-based  model  the  necessary  bond  stiffness  to  match  the  energy 
associated  with  purely  deviatoric  deformation.  The  two  can  be  combined  for  eq.  (4.13),  a  force 
state  that  clearly  indicates  the  separate  responses  to  dilation  and  deviatoric  deformation: 


t  = 


m 


m 


(4.13) 


A  quick  examination  shows  that,  in  the  case  that  the  dilation  O^d  is  not  constant,  the  force  at 
either  end  of  a  bond  will  not  satisfy  eq.  (4.5).  Thus  such  a  model  is  not  possible  in  a  bond-based 
framework  without  significant  modification. 

4.4.3  Correspondence  Models 

One  way  to  create  a  peridynamic  material  model  is  to  start  from  a  material  model  in  classical 
dynamics.  Classical  models  based  on  the  deformation  gradient  have  the  advantage  of  decades 
of  development  and  tens  of  thousands  of  hours  of  testing,  and  they  enjoy  widespread  use  in  the 
continuum  mechanics  community.  Peridynamic  correspondence  models  use  the  relative  positions 
of  a  points  neighbors  to  determine  F(Y),  a  nonlocal  approximation  of  the  deformation  gradient  F. 


F(Y) 


[  msi)(Y{«)®«)  dVf 

Jn 


K 


-1 


(4.14) 
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with  the  shape  tensor  K  defined  by 


K  =  /  u;m{^0^)dV^. 

Ju 

When  F  is  constant,  F(Y)  is  exactly  equal  to  F.  If  the  classical  model  in  question  is  hyperelastic 
with  energy  density  f2(F),  it  is  a  simple  matter  to  force  the  peridynamic  model  to  have  identical 
energy  by  defining 

W^(Y)  =  n(F(Y)) ,  (4.15) 


and  find  the  force  vector  state  by  taking  the  Frechet  derivative  of  W  with  respect  to  Y.  Alternately, 
the  classical  continuum  model  can  be  applied  to  find  the  first  Piola-Kirchhoff  stress  P  associated 
with  F: 


af2(F) 

df 


(4.16) 


The  resulting  force  state  is  calculated  from  the  stress  according  to 


T{«)  =  . 


(4.17) 


For  homogenous  deformations,  the  result  is  a  peridynamic  model  that  exactly  reproduces  the  clas¬ 
sical  model  without  ever  taking  a  derivative.  For  problems  with  very  inhomogenous  (on  the  scale  of 
the  peridynamic  horizon)  deformations,  the  peridynamic  model  will  exhibit  scale  effects  not  seen  in 
the  classical  model,  acting  to  smooth  out  the  effect  of  short-scale  deformations.  For  discontinuous 
deformations,  the  classical  model  cannot  be  evaluated  at  all,  but  the  peridynamic  correspondence 
model  will  have  no  such  problem.  It  may  be  necessary  however  to  revisit  the  choice  of  model  or 
implement  some  damage  condition. 

All  three  of  these  models  are  based  on  solid  materials.  In  such  materials,  the  fundamental 
deformation  mode  is  stretch  or  extension.  Consider  instead  a  thin  beam  deflecting  under  transverse 
load;  at  the  scale  of  the  whole  beam,  the  deflection  behavior  is  the  result  of  bending  deformation. 
To  model  the  beam  with  a  solid  material  model,  it  is  necessary  to  use  a  much  smaller  scale  -  one 
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in  which  the  beam  can  be  seen  to  stretch  on  one  side  and  compress  on  the  other.  Alternatively,  in 
a  model  whose  fundamental  mode  of  deformation  is  bending,  the  same  beam  could  be  modeled  at 
the  same  scale  as  its  behavior. 
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Chapter  5:  MODEL  DEVELOPMENT 


5.1  Bond  Pair  Material  Model 

Consider  the  material  model  illustrated  in  fig.  5.1  in  which  every  bond-vector  originating  from  a 
point  is  connected  by  a  rotational  spring  to  its  opposite  originating  from  that  same  point.  If  we  call 

. und’efdfme'd'bohd  pair . 

Figure  5.1:  Illustration  of  a  bond  pair  model  that  resists  angular  deformation 

the  deformed  angle  between  these  bonds  9,  and  choose  the  potential  energy  of  that  spring  to  be 
w{^)  =  a;(^)Q;[l  -|-  cos(6*)]  for  the  bond  pair  ^  and  — we  can  recover  the  non-ordinary  force  state 
proposed  by  Silling  in  [66]  by  taking  the  Frechet  derivative.  For  the  derivation  and  a  description 
of  the  Frechet  derivative  see  appendix  A. 


T(^)  =  Vu;(Y(^))  , 


-«  Y(^) 

|Y(^)l|Y(OI 


X 


[Y(0  Y(-^)l 

L|y(OI  |y(-^)U  ■ 


(5.1) 


Though  it  looks  complex,  eq.  (5.1)  indicates  a  bond  force  perpendicular  to  the  deformed  bond  and 
in  the  plane  containing  both  the  deformed  bond  and  its  partner  as  illustrated  in  fig.  5.2.  The  force 
magnitude  is  proportional  to  the  sine  of  the  angle  between  the  bonds  divided  by  the  length  of  the 
deformed  bond.  This  response  is  consistent  with  the  idea  of  a  rotational  spring  between  bonds  as 


Figure  5.2:  Deformation  and  force  vector  states 


long  as  the  change  in  angle  is  small.  Because  the  potential  energy  and  force  states  are  functions 
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of  pairs  of  peridynamic  bonds,  we  will  call  this  formulation  a  bond-pair  model.  Other  choices 
for  the  bond-pair  potential  function,  such  as  w  =  (tt  —  6*)^,  are  also  possible,  but  result  in  more 
mathematically  complex  analysis. 

5.2  Bond  Pair  Beam  in  Bending 

The  simplest  application  of  our  bond-pair  based  peridynamic  model  is  that  of  fig.  5.3,  a  beam  in 
transverse  bending.  Much  of  the  material  in  this  section  can  also  be  found  in  [54]. 


2(5 


Figure  5.3:  A  continuous  peridynamic  beam  with  horizon  5 


5.2.1  Energy  Equivalence 


To  determine  an  appropriate  choice  of  a  for  eq.  (5.1),  we  desire  our  peridynamic  model  to  have  an 
equivalent  strain  energy  density  to  a  classical  Euler-Bernoulli  beam  in  the  local  limit,  i.e.  when  the 
nonlocal  length  scale  vanishes.  We  will  begin  with  the  assumptions  from  Euler  beam  theory:  the 
length  of  the  beam  is  much  greater  than  thickness,  vertical  displacements  are  small,  and  rotations 
are  small.  For  small  vertical  displacements  (i.e.  sin  0  ~  ^)  we  have 


v(x  +  () 


2v(x)  -h  v(x 
~1 
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where  v  is  the  vertical  displacement  of  material  point.  Momentarily  assuming  that  v  is  continuous 
and  using  a  Taylor  series  to  expand  the  right-hand- side  of  eq.  (5.2) 


^(Y(0,Y(-0)^7r-e^  +  O(a, 

TT  -  +  0{^^)  , 


(5.3) 


dx^  ■ 

Substituting  eq.  (5.3)  into  the  equation  for  the  strain  energy  density  of  a  single  bond-pair, 

u;(0  =a;(e)«[l  +  cos(^(Y(0,Y(-0))]  , 

If  we  use  a  weighting  function  C(j(^)  =  a;(|^|)  and  assume  that  the  u  plays  the  role  of  a  localization 
kernel,  i.e.  cu  =  0  V  ^  >  <5,  the  resulting  strain  energy  density,  W,  for  any  material  point  in  the 
peridynamic  beam  is 

Equating  W  with  the  classical  Euler-Bernoulli  beam  strain-energy  density,  D,  and  taking  the  limit 
as  (5  — )■  0  we  can  solve  for  a 


lim  IE  =  , 

(5^-0 

a  2  2 

-ruK  =  —K  , 

El 

a  =  —  , 
m 


(5.4) 


with 


m=  a;(^)C^d^. 

J-s 
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While  this  demonstrates  the  model’s  equivalence  to  a  linearly-elastic  Euler  beam,  if  we  keep 
an  additional  term  from  the  Taylor  series  approximation  of  eq.  (5.2),  we  recover  a  slightly  more 
complex  expressions  for  change  in  angle  that  is  demonstrated  in  5.2.2  to  reproduce  an  Euler  beam 
governed  by  Eringen’s  model  of  nonlocal  elasticity. 

5.2.2  Relation  to  Eringen  Nonlocality 


If  we  keep  an  additional  term  from  the  Taylor  series  approximation  of  eq.  (5.2),  we  recover  a 
slightly  more  complex  expressions  for  change  in  angle 


0mO,Y{-C))  ^  arctan  (  tt  - 


and  for  the  strain  energy  (again  substituting  n  =  v”  for  readability). 


As  the  horizon  5  becomes  small,  higher-order  ^  terms  become  relatively  less  important,  and 
is  dominated  by  for  large  k  and  by  ^'^kk"  for  small  k.  The  remaining  terms  can  be  rearranged. 


Z  D 


'-s 


in  a  manner  strongly  suggesting  an  alternative  bending  resistance  term.  We  can  picture  a  bending 
resistance  based  on  the  bond  length  and  proportional  to  the  nonlocal  curvature  k  =  (k  +  so 

that 


K  = 

w  ^ 


r-S  ^2 

/  u)(^)a—KKd^. 

J-s  2 


(5.5) 
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The  same  analysis  can  be  taken  further  to  obtain  higher-order  energy  terms  with  even  powers  of 
^  and  even  order  derivatives  of  k.  Not  all  of  these  higher-order  terms  can  be  separated  into  the 
product  of  a  local  curvature  and  nonlocal  bending  resistance. 

Eringen’s  model  for  nonlocal  elasticity  in  [21]  begins  with  a  nonlocal  modulus  (denoted  here 
as  K{\x.'  —  x|,  r))  that  relates  the  nonlocal  stress  t  at  a  point  to  the  classical  (local)  stress  cr  in  the 
nearby  material  through  the  integral 

t=  f  /f(|x' —  x|,  r)(T(x')du(x') . 
iv 

In  the  local  limit  these  relationships  take  the  form  of  higher-order  gradients.  Using  a  1 -dimensional 
decaying  exponential  nonlocal  modulus  K{\x\,t)  =  results  in  a  relationship  between  fiD 

and  (Tid, 


in  which  (rl)  is  a  scale-based  material  parameter.  For  well-behaved  fiD  and  ctid  and  small  values 
of  a'jD  and  (t()^,  we  can  see  that  this  relationship  could  be  reformulated  as 

fiD=  +  aiD. 

If  we  consider  the  results  of  the  previous  section  and  let  dM  =  yadA  and  a  —  Eyn,  the  contribu¬ 
tion  to  moment  resulting  from  Eringen’s  nonlocal  elasticity  in  a  fiber  at  y 

Ey‘^{K  +  {tIYk")  ,  (5.6) 


and  the  resulting  strain  energy 


b{y)E^K{K  +  {TlfK'')dy , 
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bear  a  striking  resemblance  to  eq.  (5.5).  In  fact,  by  carefully  choosing  peridynamic  parameter 
values,  the  results  can  be  made  identical.  For  a  rectangular  beam  of  width  b  and  thickness  t, 
choosing 


S  =  , 


Ebt^ 
54(t()^  ’ 


results  in 


the  same  result  for  both  models. 

The  similarity  between  eqs.  (5.5)  and  (5.6)  is  not  accidental;  Eringen’s  gradient  elasticity  is 
the  solution  to  the  integral  formulation  of  the  nonlocal  stress  integral  equation  just  as  the  peridy¬ 
namic  energy  is  an  integral  function  of  nonlocal  displacements.  It  is  therefore  unsurprising  that, 
like  Eringen’s  nonlocal  elasticity  [12],  this  peridynamic  bending  model  fails  to  predict  the  stiffen¬ 
ing  associated  with  nanoscale  cantilevers.  Instead,  the  advantage  of  peridynamic  models  is  their 
natural  handling  of  discontinuities. 


5.2.3  Weighting  function  and  inelasticity 

The  weighting  function  ix){^)  describes  the  relative  contribution  of  each  bond-pair,  and  can  be  de¬ 
fined  according  to  physical  or  mathematical  considerations.  While  any  function  that  produces 
a  convergent  integral  for  m  will  reproduce  an  elastic  Euler  beam,  a  physically  meaningful  choice 
of  Lo  will  allow  us  to  extend  our  model  to  certain  inelastic  behaviors.  Consider  a  classical  Euler- 
Bernoulli  beam  in  bending  with  curvature  k.  Fibers  running  parallel  to  the  neutral  axis  of  the  beam 
are  stretched  in  proportion  to  their  distance  from  the  neutral  axis,  with  strain  e  =  yn.  If  the  fibers 
are  linearly  elastic,  then  the  axial  stress  at  each  location  is  cr  =  =  Eyn,  and  the  contribution 

to  supported  moment  dM  =  nEy'^dA.  By  comparing  the  formulations  for  the  moments  carried 
by  the  Euler  beam  in  fig.  5.4  and  those  of  the  bond-pair  beam  in  fig.  5.5,  we  see  some  definite 
parallels. 
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t 


y  =  V2 


y  =  -V2 

a  =  Ee 


dM  =  yadA 


Figure  5.4:  Euler  beam  moment  contribution 


-5 


p'—ij^  6  ^  '  dM  ^  sm{2d)u>{^)d^ 

Figure  5.5:  Bond-pair  moment  contribution 


dM  =  ^ed^ 


Me  = 
MpD  = 


a  y  dA 


l-s 

sin(A^) 


Ek  b{y)dy , 


(5.7a) 


i-s 


(5.7b) 


The  term  y  is  the  distance  from  the  beam’s  neutral  axis  and  h{y)  is  the  width  of  the  beam  at 
that  distance  from  the  neutral  axis.  The  similarity  between  classical  and  peridynamic  moment 
formulations  in  eqs.  (5.7a)  and  (5.7b)  suggests  a  possible  formulation  for  the  weighting  function: 


at  y  =  ^\- 


(5.8) 


This  weight  function  analogizes  the  relative  contributions  of  bond  pairs  of  different  lengths 
to  the  relative  contributions  of  fibers  at  different  distances  from  the  centerline.  An  example  for  a 
rectangular  beam  is  illustrated  in  fig.  5.6.  For  an  I  beam  with  height  /ibeam>  width  rcbeam>  web  height 
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Figure  5.6:  Weight  function  for  a  beam  of  rectangular  cross-section 


/iweb>  and  web  width  w^sb,  substituting  the  beam  profile 


'^web  if  I?/ 1  — 


h 


web 


2  ’ 


Kv)  =  { 


.  r  ^web  ^11^  ^beam 

-(^beam  it  <  \y\  < 


0 


otherwise , 


into  eq.  (5.8)  gives  the  weight  function 


^(0 


I^IWweb 

^beam 

^  l^kbeam  if()^^  <  |^|  <  <)  , 

^beam 

0  otherwise , 


and  is  ilustrated  in  fig.  5.7.  While  this  weighting  function  offers  no  advantages  over  a  uniform 
weight  function  in  the  case  of  the  linearly  elastic  beam,  it  offers  a  way  to  model  advancing  plastic¬ 
ity. 

In  a  deformed  elastic  perfectly-plastic  beam,  axial  fibers  are  still  stretched  in  proportion  to  their 
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Figure  5.7:  Weight  function  for  an  I-beam 


distance  from  the  neutral  axis,  but  the  relationship  a  =  Ee  =  Eyn  only  holds  for  |e|  =  \yK\  <  Cc- 
For  greater  stretches,  the  relationship  becomes  a  —  ±Eec.  To  model  this  behavior,  consider  a 
bond  pair  with  similar  behavior:  for  angular  deformation  less  than  some  critical  angle,  the  model 
behaves  as  previously  described,  but  the  magnitude  of  the  force  remains  constant  above  a  critical 
deformation  according  to 


|T(OI 


sm(^(X(0.X(-0)) 

|Y(OI 


if^  <  Oc, 


au{^) 


sin(6c) 

|Y(OI 


if  0  >  6c. 


(5.9) 


To  determine  the  critical  angle  dc,  we  let  the  onset  of  plasticity  in  pairs  of  the  longest  bonds  to 
coincide  with  the  onset  of  plasticity  in  the  fibers  at  the  top  and  bottom  surfaces  of  the  classical 
beam.  For  small  curvatures  AO  =  A6c  =  For  curvatures  |«:|  >  Kc  =  the  radius 

within  which  bonds  are  in  the  elastic  region  is  (ie  =  and  parallels  the  distance  from  the  beam 
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centerline  that  fibers  are  in  the  elastic  region  y, 


t  ^ 

2  K 


rve 


^classical  =  2  /  Eb{y)y  Kdy  +  2  /  Eb{y)ecydy  , 

Jo  J  ye 

MpD  =  2  f  a(jj{^)^‘^Kd^  +  2  f  au{^)A0cCdC  ■ 
Jo  JSe 


Of  course,  as  long  as  the  force  is  independent  of  history,  this  model  only  represents  a  nonlinear 
elastic  material.  By  keeping  traek  of  the  plastic  deformation  =  9  —  OcOf  each  bond-pair,  and 
applying  it  as  an  offset,  we  can  reproduce  the  hysteresis  associated  with  elastic-perfectly-plastic 
deformation. 

More  simply,  we  can  model  a  brittle  material  by  setting  the  force  to  zero  for  bond  pairs  exeeed- 
ing  a  critical  angle. 


mo\ 


sin(e(Y(0,Y(-0)) 

|Y(?>I 


if  ^  ) 


(5.10) 


0 


if  9  >  9c, 


and  additionally  recording  bond  pairs  that  have  exceeded  their  eritical  angle  and  permanently  set¬ 
ting  their  influence,  i.e.  co,  to  zero. 


5.3  Bond  Pair  Plate  in  Bending 

The  next  case  we  will  analyze  is  the  extension  of  the  bond  pair  beam  model  to  fig.  5.8,  a  flat  plate 
in  the  xy  plane,  with  displacement  w  in  the  ;2-direction.  Much  of  the  material  in  this  section  can 
also  be  found  in  [55]. 

5.3.1  Energy  Equivalence 

As  with  the  beam  model,  we  determine  an  appropriate  choiee  of  a  so  that  our  peridynamic  model 
will  have  an  equivalent  strain  energy  density  to  a  classical  Kirckhoff  plate  in  the  local  limit.  We 
will  begin  with  the  assumptions  from  Kirckhoff  plate  theory:  straight  lines  normal  to  the  mid¬ 
surface  remain  both  straight  and  normal  to  the  deformed  mid-surface,  and  the  plate  thickness  does 
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Figure  5.8:  Illustration  of  a  bond  pair  on  a  plate. 


not  change  with  deformation.  As  with  the  Euler  beam  energy  equivalence,  we  will  start  with  the 
original  assumptions  from  Kirchhoff-Love  plate  theory  of  small  displacements  and  rotations,  but 
they  will  not  constrain  the  validity  of  the  model  for  larger  displacements  and  rotations.  For  small 
vertical  displacements  we  have 

where  w  is  the  vertical  displacement  of  material  point.  Taking  ^  =  ^(cos(0),  sin((/)))  in  cartesian 
coordinates  and  momentarily  assuming  continuous  displacements  for  the  sake  of  comparison,  we 
use  a  Taylor  series  to  expand  the  right-hand- side  of  eq.  (5.1 1)  about  C  =  0 


Y(-^))  I  (cos^(0)«:i  -b  sin^(0)fi;2  +  2  sin(^)  cos(0)«:3)  -b  ,  (5.12) 


with 


Ki 


d^w  d^w 

dx\ ’  dxl 


d‘^w 

dxidx2 
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Substituting  eq.  (5.12)  into  the  equation  for  the  strain  energy  density  of  a  single  bond-pair, 


u;  =  a;(O«[l  +  cos(0(Y(^),Y(-O))]  , 

=  a;(^)Q;— cos^(0)  -|-  sin^(^)  -|-  2k,iK,2  cos^(^)  sin^(^)  -|-  cos^(^)  sin^(^) 

8 

+  4kiK3  cos^(^)  sin(^)  -|-  4:K2K3  cos(0)  sin^(0))  -|-  . 

If  we  use  a  weighting  function  a;(^)  =  co{^)  and  assume  that  the  co  plays  the  role  of  a  localization 
kernel,  i.e.  a;  =  0  V  ^  >  5,  the  resulting  strain  energy  density,  W,  for  any  material  point  in  the 
peridynamic  plate  is 

n27r 

w^d^d^, 

_ 

=«—  +  ^2  +  -KiKi  +  2^3^  J  +  0{5^)  . 

Equating  W  with  the  classical  Kirchhoff  plate  strain-energy  density,  f2,  and  taking  the  limit  as 
(5  — )■  0  we  can  solve  for  a 


lim  kk  =  O  , 
(5^-0 

StT  /^2,2,2  _L^2\ 

a—m  Ik^  +  K2  +  -KiKi  +  -K^j  = 


jih^ 


a  = 


[12(1-1.) 

2/r/i3 

3m  ’ 


(ki  +  k\  +  2vk.iKi  +  2(1  —  i')k^ 


J  u=\l?, 

(5.13) 


with 


n27r 


where  fi  is  the  shear  modulus,  h  is  the  thickness  of  the  plate,  and  we  have  evaluated  the  classical 
Kirchhoff  strain-energy  at  a  Poisson  ratio  of  ^/s  in  order  to  solve  for  alpha  as  a  constant.  Because 
a  is  inversely  proportional  to  m,  the  energy  does  not  change  with  varying  choices  for  lo  and  5. 
It  should  be  noted  that  the  restriction  v  =  i/s  is  the  same  imposed  by  the  use  of  a  bond  based 
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peridynamic  model  for  in-plane  deformation  of  a  2D  peridynamic  plate.  We  will  show  an  extension 
to  this  model  that  removes  this  restriction  in  Section  5.4. 

5.3.2  Combining  Bending  and  Extension  Models 

The  bond-pair  bending  model  does  not  resist  in-plane  stretching  or  shear  deformation  because 
these  deformations  preserve  the  angles  between  opposite  bonds.  If  these  behaviors  are  expected  in 
combination  with  bending,  a  useful  model  must  resist  both  in-plane  and  transverse  deformations. 
To  create  a  plate  model  that  also  resists  these  deformations,  i.e.  a  flat  shell,  we  combine  the  bond- 
pair  model  with  a  two-dimensional  version  of  the  original  bond-based  linearly-elastic  peridynamic 
solid  model  from  [64].  In  this  model,  individual  bonds  act  as  springs  resisting  changes  in  length; 

T«>=/J(|Y{?)|-|«|)i^.  (5.14) 

By  matching  the  energy  of  a  2D  material  in  shear  deformation,  we  can  relate  (3  to  the  shear  modulus 
and  thickness  of  the  shell.  Following  the  example  of  [66],  we  begin  with  a  2D  material  under  pure 
in-plane  shear.  In  Einstein  notation,  the  strain  energy  of  this  material  is 

ICpD  = 


where  is  the  deviatoric  strain  tensor.  To  evaluate  the  integral  we  will  exploit  the  symmetry 
properties.  With  i,j,  fc,  /  =  1, 2.  For  a  circular  uj{^)  =  c<;(|^|),  combinations  of  {i,  k,  /}  with  an 
odd  number  of  each  index,  such  as  {1,1, 1,2}  or  {2, 1,2, 2),  will  result  in  odd  powers  of  sine  and 


^  h  , 


dA^  , 


2  7a 
/  ‘^(0 


•2Ja'"'  l«l  l«l 
^44  {  d.4s  , 


I A 
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cosine  and  integrate  to  0. 


m  =  I 

[3(eiieii  +  622^22)  +  (^11622  +  ei2ei2  +  ei2e2i  +  £21^12  +  ^21621  +  622^11)] , 


VVpY,  — 


'A 

P  m  I 

fd  d  d  /c  X  XX  X  X  ^ 

^ij^klK^ij^kl  +  OjfcOji  +  Oiidjk)  , 


P  m 


e'J'e'J*. 

'-ij  '-ij 


P  = 


8  fi  h 


m 


Having  calibrated  the  bond-extension  model  to  the  shear  modulus  for  a  case  of  pure  in-plane  shear, 
applying  a  different  uniform  strain  (such  as  might  result  from  uniaxial  tension)  reveals  the  bond- 
based  model  to  result  in  a  one-parameter  linearly-elastic  model  with  Poisson’s  ratio  u  =  1/3. 

Combining  the  bending  and  extension  models  allows  for  the  description  of  more  complex  be¬ 
haviors,  particularly  the  stiffening  effect  of  in-plane  tension  on  the  transverse  bending  of  a  shell. 
Consider  a  single  bond-pair  in  the  combined  model  shown  in  Fig.  5.9.  As  the  two  sides  are  pulled 


^')  extension 
—  bending .  — 


T(0+T(-0 


-T(0 


extension 


bending 


Figure  5.9:  The  Hybrid  Model  Combines  Bending  and  Extension  Components 

apart,  the  magnitude  of  the  extension  force  in  each  bond  increases,  and  the  magnitude  of  the  bend¬ 
ing  force  decreases.  At  the  same  time,  the  angle  at  which  the  extension  force  acts  decreases,  and 
the  angle  of  action  for  the  bending  force  increases.  For  small  amounts  of  bending  and  reasonable 
stretches,  increased  tension  in  the  direction  of  the  bond  pair  results  in  increased  restorative  force. 
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5.4  Extension  to  arbitrary  Poisson  ratio 


Although  many  materials  have  Poisson  ratios  of  i/  ~  ^/s,  it  is  nonetheless  desirable  to  extend  the 
model  to  materials  with  arbitrary  Poisson  ratios.  For  isotropic,  linearly  elastic  models  of  solid 
materials,  Silling  et  al.  extended  the  peridynamic  material  model  to  arbitrary  material  parameters 
in  [66]  by  decomposing  the  deformation  into  isotropic  and  deviatoric  components.  In  the  absence 
of  plastic  deformation,  we  need  only  find  the  difference  between  the  strain  energy  of  a  deformed 
bond-based  plate  and  the  strain  energy  of  an  elastic  plate  with  Poisson’s  ratio  v  ^  1/3.  The  differ¬ 
ence  is  a  function  of  the  isotropic  strain  in  two  dimensions, 


W*  = 


fj,  h  f  Sp  —  1 


^2D  1 


O2D  = 


^^total 


m  Ja 
ji  h  f  3p 


i-p 

^m\mo\-\^\)dA^, 


1 


i-p 


^20  + 


4:  fX  h 


m 


This  is  to  be  expected  because  the  bond-based  model  was  calibrated  to  the  shear  strain  energy, 
leaving  discrepancies  proportional  to  the  isotropic  strain  energy  that  fall  to  0  as  Poisson’s  ratio 
approaches  p  =  '^f 3. 

This  decomposition  method  inspires  a  similar  approach  to  our  plate  model.  To  perform  the 
same  extension  for  the  plate  model  in  bending,  we  find  the  error  in  the  1 -parameter  strain  energy 
for  p  7^  1/3 


W* 


w* 


fib? 


12(1  -  p) 
fih^ 


+  K.I  +  2pKiK2  +  2(1  —  p)k.I) 


992  1 

\  +  ^2  +  +  2(1  —  - 


=2/r 


12(1-1) 

3p  —  1  f  Ki  +  ^ 

12  1  -  1/ 


The  discrepancy  in  energy  is  proportional  to  the  square  of  average  curvature,  which  we  will 

also  refer  to  as  the  isotropic  curvature.  The  isotropic  curvature  can  be  envisioned  as  the  portion  of 
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the  deformation  that  resembles  a  hemispherical  bowl.  The  remainder  of  the  bending  deformation, 
that  which  is  left  when  the  isotropic  curvature  is  subtracted  out,  resembles  a  saddle.  This  remaining 
component  is  the  deviatoric  deformation,  and  both  components  are  shown  in  fig.  5.10.  Note  that 
the  orientation  of  the  deviatoric  bending  will  change  depending  on  the  particular  curvature  being 
decomposed,  while  the  isotropic  curvature  will  only  change  in  scale.  A  complete  decomposition  of 


Total  Bending 


Bending  Decomposition 
Isotropic 


Deviatoric 


bending  energy  into  isotropic  and  deviatoric  components  as  performed  by  Fischer  in  [24]  produces 
a  far  more  complex  model  and  is  unnecessary  at  this  time.  For  a  single  bond  pair  we  can  represent 
the  curvature  vector  along  the  bond  pair  as 

^  -  icp 

For  large  rotations,  we  can  define  an  average  curvature  vector  R,.  This  leads  us  to  model  the  average 
curvature  as 


1 


S  p2it 


"“mX  y„ 

n27r 


Y(0+Y(-0 


ed(/>de, 


The  weighting  function  performs  the  same  function  as  in  the  previous  section.  We  can  rewrite 
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the  energy  discrepancy  in  terms  of  n. 


TTT*  ^  3p  —  I  _o 

W  =  2/1—- - 


12  1 


We  can  take  the  Frechet  derivative  (details  in  A)  to  produce  a  correction  force  vector  state 


T!«) 


8/U  3iy  —  1  uj{^)  _ 
m  12  1  —  z/ 


(5.15) 


that  is  not  directly  dependent  on  the  deformation  of  a  single  bond  pair.  Instead,  eq.  (5.15)  represents 
a  bond-length  dependent  “pressure”  applied  to  every  pair  of  bonds  extending  from  a  node.  This 
“pressure”  is  proportional  to  the  curvature  vector  at  that  node.  A  weighting  function  co{^)  =  |^| 
can  ensure  that  the  integral  expression  for  force  at  a  point  is  convergent.  This  extra  term  that 
is  dependent  on  the  bending  of  all  the  pairs  around  a  material  point  means  that  the  extension  is 
not  properly  a  bond-pair  model.  Instead,  it  would  be  more  accurate  to  call  it  a  bond-multiple 
model,  in  which  the  bond  forces  and  energies  are  functions  of  the  relationship  between  a  family  of 
bonds.  In  either  the  continuous  or  discrete  cases,  this  model  extension  requires  the  additional  step 
of  evaluating  the  isotropic  curvature  at  each  point,  but  the  increased  complexity  of  the  extended 
model  captures  in  the  local  limit  the  behavior  of  a  two-parameter  elastic  material  plate. 
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Chapter  6:  NUMERICAL  SIMULATION 


The  models  developed  in  chapter  5  are  all  based  on  a  continuum  material  containing  an  infinite 
number  of  points.  The  model  parameters  and  behaviors  are  defined  by  integrals  over  length  or 
area,  and  only  make  sense  if  properties  are  defined  continuously.  This  is  not  to  say  that  the  prop¬ 
erties  themselves  need  be  continuous,  but  they  cannot  be  defined  only  at  a  finite  number  of  points. 
While  such  models  are  mathematically  convenient,  they  are  less  convenient  for  performing  many 
engineering  problems.  For  the  analysis  of  these  problems,  it  is  useful  to  have  a  discretized  version 
of  the  model  that  can  be  evaluated  from  the  properties  of  a  finite  number  of  peridynamic  nodes  that 
represent  the  continuous  feature. 


Figure  6.1:  Translating  Continuum  Peridynamics  to  a  Discrete  Domain 


Discretizing  the  continuous  domain  of  interest  as  illustrated  in  fig.  6.1  is  a  straightforward 
task.  First,  the  entire  domain  is  divided  into  small,  non-overlapping  subdomains.  Ideally,  these 
subdomains  are  regular  in  size  and  compact  in  shape,  and  it  is  important  that  the  largest  dimen¬ 
sion  of  a  subdomain  be  similar  to  the  smallest  dimension  as  well  as  significantly  smaller  than  the 
peridynamic  horizon. 

A  peridynamic  node  is  then  placed  at  the  centroid  of  each  subdomain.  This  node  has  no  volume 
(is  a  zero-dimensional  point),  but  it  does  have  mass  equal  to  the  mass  of  the  material  in  the  subdo¬ 
main.  After  the  nodes  are  placed,  neighborhoods  are  constructed.  The  neighborhood  of  the  node 
at  Xj  consists  of  all  the  nodes  Xj  for  whom  |xj  —  Xj|  is  smaller  than  the  horizon.  As  the  number 
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of  nodes  in  each  neighborhood  increases,  the  material  represented  by  the  nodes  in  a  neighborhood 
will  converge  to  the  material  within  the  horizon  in  the  continuous  domain.  At  the  same  time,  the 
integrals  of  well-behaved  functions  over  neighborhood  regions  in  the  continuous  model  can  be 
approximated  by  weighted  sums  of  the  same  functions  evaluated  at  the  discrete  nodes  in  a  neigh¬ 
borhood.  This  is  the  key  to  translating  the  continuous  peridynamic  model  into  a  more  practical 
discrete  version. 

6.1  Discretized  Bond  Pair  Beam 

Discretizing  the  bond-pair  model  is  primarily  matter  of  exchanging  integrals  for  sums, 

wi$,)=u{^,)a[l  +  cos{e{Y{^,),Y{-m]  , 

_  />  f  ^(x  +  -  2t;(x)  +  t;(x  - 

I?  J  2  V 

in  which  is  the  z*  bond  emanating  from  the  point  x  to  each  of  the  n  points  within  distance  5  of 
point  X. 


c  Ar 

a  = - ;  c  =  El;  m  = 

m  ^ 

w  =  Ax  j2  , 

i=\  ^  V  / 

Discretization  of  the  original  model  results  in  the  equation  of  motion. 


p(x)u(x)  =  f(x)  + 


«(x)  Pi 


X 


p*l  IP* 

«(x  +  ^i)  (-Pi) 


Pi  I  |qi| 


X 


Ipil 


-Pij 


X 


Ipil 


(6.1) 


(6.2) 
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with 


Pi  =  $i  +  u(x  +  -  u(x) , 

q*  =  -$i  +  u(x  -  $i)  -  u(x) , 
ri  =  $i  +  u(x  +  2^i)  -  u(x  +  , 

and  for  small  displacements  and  rotations  in  a  uniform  beam, 

p(x)v(x)  =/(x) 

+«  2a;(f )  ( ~  ~  ~  ^  ~  ^  ^  ^ 

i  ^ 

It  is  worth  noting  the  similarity  between  this  expression  and  a  finite-difference  fourth  derivative  of 
displacement,  a  result  expected  from  Euler  beam  theory.  This  discretization  requires  that  nodes  be 
evenly  spaced  along  the  entire  beam,  otherwise  the  displacement  i>(x  —  is  ill-defined.  For  this 
reason,  the  discretization  does  not  allow  for  areas  of  higher  and  lower  “resolution”. 

6.2  Discretized  Bond  Pair  Plate 

As  with  the  beam,  discretizing  the  bond-pair  model  is  primarily  matter  of  exchanging  integrals  for 
sums; 


c  (Ax)^ 

a  = - ;  c  = 


m  = 


i=l 


)r 

2(1  -  i')  12  V  ifil 


(l-j/)12 

W  =  ( io(x  +  ^.)-2^(x)+^x-i.) 


(6.3) 


Discretization  of  the  1 -parameter  bending  model  results  in  the  same  equation  of  motion  as  for  the 
beam  model  (eq.  6.2). 
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Figure  6.2:  Discretized  peridynamic  plate  with  illustrated  bond  pair 


Implementing  the  2-parameter  model  requires  finding  the  isotropic  curvature  at  each  point. 


-  /  X  1  /'/•  nP*  + 

.2  > 

m  F 

Sz 


/‘“(x)  =  ^  I  [a“(x)ft(x)  -  a‘”(x  +  ^j)«(x  +  ^j)] 


£2 


As  with  the  discretized  beam,  the  discretization  of  the  bond-pair  plate  (fig.  6.2)  must  be  absolutely 
regular.  Discretizing  the  bond-pair  model  as  proposed  above  requires  that  nodes  be  evenly  spaced, 
Ax,  throughout  the  entire  plate,  otherwise  the  displacement  w{'x  —  £j)  is  undefined.  For  this 
reason,  the  discretization  does  not  allow  for  areas  of  higher  and  lower  “resolution”.  This  restriction, 
while  inconvenient  in  the  ID  case,  is  fairly  restricting  for  plate  analysis.  An  extension  to  this 
discretization  that  would  allow  changing  mesh  resolution  will  require  interpolation  between  the 
nodes. 
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6.3  Numerical  Model  Extensions 


6.3.1  Curved  Shapes 

On  a  curved  surface,  the  location  of  the  point  x  —  ^  might  be  off  of  the  surface  entirely.  One 
method  of  applying  the  bond-pair  model  to  eurved  surfaces  is  through  the  use  of  “virtual”  points. 
These  points  have  no  mass  and  do  not  have  families  of  peridynamic  neighbors,  they  only  allow  the 
definition  of  bond  pairs  that  are  straight  in  the  undeformed  eonfiguration.  In  the  simplest  method. 


Figure  6.3:  Virtual  Points  Allow  Straight  Pairs  on  Curved  Surfaces 

each  virtual  point  is  located  just  above  or  below  a  real  point  in  the  model.  In  this  case,  properties 
sueh  as  displaeement  are  taken  to  be  the  same  as  for  the  nearby  real  point.  Beeause  the  virtual 


x(-«)  . . , 


Figure  6.4:  Virtual  Points  Take  the  Displacement  of  Nearby  Real  Points 

point  has  no  mass  is  not  part  of  any  other  bond  pairs,  it  cannot  be  assigned  a  foree.  Instead,  the 
foree  on  a  virtual  point  resulting  from  deformation  of  a  bond  pair  is  instead  applied  to  the  nearest 
real  points.  This  results  in  a  straightforward  extension  of  the  bending  model  from  flat  plates  (and 
beams)  to  features  that  have  curvatures  that  are  small  over  the  peridynamic  horizon. 
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6.3.2  Irregular  Discretization 


A  curved  surface  is  not  the  only  reason  to  implement  virtual  points,  and  even  many  curved  surfaces 
do  not  allow  for  regular  discretization.  When  discretization  is  irregular,  due  to  three-dimensional 
curvature,  irregular  shapes,  or  a  need  for  increased  resolution  in  some  areas,  there  are  necessarily 
points  at  which  there  is  no  real  point  at  the  location  of  x  —  An  example  of  changing  mesh 
density  resulting  in  a  need  for  interpolation  can  be  found  in  fig.  6.5,  which  shows  a  small  family  of 
nodes  at  the  edge  of  a  change  in  discretization  coarseness.  Note  that,  while  bonds  p2  and  q2  form 


Figure  6.5:  Virtual  Points  Pair  Up  Unpaired  Neighbors 

a  perfect  bond  pair,  there  is  no  bond  exactly  opposite  pi.  To  solve  this,  we  add  a  virtual  point  to 
create  a  bond,  gi,  that  will  form  a  pair  with  pi.  Because  this  point  is  not  part  of  the  discretization, 
it  has  no  mass,  and  its  properties  must  be  determined  from  the  properties  of  the  surrounding  nodes. 

An  easy  method  of  determining  properties  (such  as  displacement)  at  virtual  nodes  is  to  use  a 
weighted  average.  For  an  irregular  straight  beam,  determining  the  values  of  properties  at  virtual 
points  is  simple.  To  determine  the  value  of  a  property  at  point  C,  we  used  a  weighted  average  of 
the  values  of  that  property  at  the  nearest  two  real  points,  A  and  B.  The  weight  value  of  5  is 
determined  to  make  a  linear  interpolation  (or  extrapolation)  using  the  ^-coordinates  of  the  nearest 


58 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


two  real  points. 


Wb  = 


a 


~  Bx 


1 


Wa  —  1  —  wb  • 


The  problem  becomes  a  little  more  complicated  for  curved  beams,  but  we  can  tackle  it  by  project¬ 
ing  the  virtual  point  onto  the  line  between  the  two  nearest  real  points.  The  weighting  function  is 
then  computed  according  to 


AC  = 


BC  = 


Wb={ 


AB 

^  AB 

■AC 

\AB\  ' 

BA 

r  BA 

■  BC 

\BA\  ' 

\Ad\ 

\AB\ 


l^g| 

I  \AB\ 

Wa  =  1-  Wb  . 


if  \BC  I  <  \AB\ , 


Otherwise , 


Determining  properties  at  virtual  points  is  more  difficult  in  plate  and  shell  models.  One  method 
of  generating  useful  weights  that  is  relatively  robust  is  barycentric  interpolation.  We  start  by  finding 
the  three  (non-colinear)  real  nodes  closest  to  the  location  of  the  virtual  node,  A,  B,  and  C.  Next,  we 
find  the  signed  areas  of  the  triangles  ABC,  ABX,  BCX,  and  CAX,  with  X  being  the  virtual  node. 
The  weight  of  node  A  is  the  area  ratio  between  BCX  and  ABC,  the  weight  of  node  B  is  the  ratio 
of  areas  CAX  and  ABC,  and  the  weight  of  node  C  is  the  ratio  of  areas  ABX  to  ABC.  Using  signed 
areas  allows  the  weights  to  be  negative  to  extrapolate  properties  of  a  virtual  node  outside  of  ABC. 
Because  these  weights  are  calculated  from  the  initial  positions  of  the  node,  they  can  be  stored  for 
swift  evaluation  of  properties  at  virtual  nodes. 

With  the  properties  of  the  virtual  points  determined,  the  model  can  be  evaluated  in  the  same 
manner  as  the  uniformly  discretized  models  of  the  previous  papers.  Where  forces  are  calculated 
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Figure  6.6:  Barycentric  interpolation  is  based  on  the  relative  areas  of  sub  triangles 

to  act  on  a  virtual  node,  those  forces  are  redistributed  to  the  supporting  real  nodes  according  to 
the  weight  each  point  has  in  the  interpolation.  Barycentric  interpolation  is  linear  and  therefore 
exactly  reproduces  the  linear  displacement  fields  in  fig.  6.7,  including  extrapolation  well  outside  the 
interpolation  points  (all  of  which  are  within  the  square  with  corners  (0,0)  and  (1,1)).  Unfortunately, 
as  fig.  6.8  demonstrates,  barycentric  interpolation  is  not  exact  for  quadratic  surfaces.  The  difference 
between  the  quadratic  surface  and  the  linear  interpolation  decreases  with  denser  discretization  as 
the  curvature  of  the  surface  between  interpolation  points  decreases,  as  demonstrated  in  fig.  6.8. 
This  method  therefore  requires  that  the  curvature  of  the  surface  be  small  relative  to  the  peridynamic 
horizon  to  ensure  accurate  virtual  node  properties. 

The  same  method  of  virtual  nodes  also  allows  the  modeling  of  curved  surfaces,  in  which  the 
perfect  opposite  of  a  bond  may  not  lie  near  but  not  on  the  surface  of  the  plate  or  shell.  As  long 
as  the  curvature  of  the  surface  is  small  (at  the  scale  of  the  peridynamic  horizon),  each  resulting 
virtual  nodes  will  be  nearly  in  the  plane  formed  by  its  nearest  neighbors.  Finding  the  weights  of 
the  surrounding  nodes  is  performed  just  as  in  the  planar  case,  except  that  the  areas  are  formed 
between  the  projection  of  the  virtual  node  location  X  onto  the  plane  formed  by  A,  B,  and  C. 

To  compute  the  weight  of  node  A  in  the  interpolation  of  properties  at  virtual  node  X,  let  AB 
represent  the  vector  from  node  A  to  node  B,  and  use 


IVa  = 


B-X 

- • 

2 


BC  X 


f  BCxBA\ 
VI^C  X  BA\) 


(6.4) 
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B  ary  centric  Interpolation  of  Plane 
Point  Estimates 


Figure  6.7:  B  ary  centric  estimate  and  error  for  plane  interpolation 


Barycentric  Interpolation  of  Quadratic  Surface 


Actual  Values 


Real  Points 


0.6  n  q^^O.2 


0  0.2  0.4  0.6  0.8  1.0 


Point  Estimates 


Error 


Figure  6.8:  Barycentric  estimate  and  error  for  quadratic  surface 
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-OTX 


Barycentric  Interpolation  of  Quadratic  Surface 

Real  Points 

Actual  Values  ^  ^  ^  I  ^ 


Figure  6.9:  Barycentric  estimate  improves  with  denser  discretization 
After  finding  Wb  and  Wc  in  similar  fashion, 


Wa  = 


Wa 

Wa  +  Wb^Wc' 


(6.5) 


If  the  projection  of  X  onto  the  plane  defined  by  A,  B,  and  C  lies  outside  the  triangle  ABC,  one  or 
two  of  W A,  Wb  and  Wc  will  be  negative,  though  they  will  still  sum  to  1. 


6.3.3  Extended  Discretization 


With  the  addition  of  virtual  points  and  the  incorporation  of  irregular  discretizations,  it  is  necessary 
to  reformulate  the  discretized  bending  model.  The  bond  pair  coefficients  at  point  x  of  a  beam  in 
bending  (eq.  (6.1))  become 

c/(x)  /(x  +  ^J  ^  c  =  EI,  m  =  ^w(^J^,2/(x  +  ^J  , 
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with  /(x)  and  /(x  +  representing  the  lengths  of  beam  represented  by  the  nodes  at  (x)  and 
(x  +  ^j).  Similarly,  the  plate  eoeffieients  from  eq.  (6.3)  beeome 


c  A(x)A(x  +  n 

= - >  c  =  - - -  ,  m 

m  (1  —  12 


i=l 


with  A(x)  and  A(x  +  representing  the  areas  represented  by  the  nodes  at  (x)  and  (x  +  ^j).  Note 
that  for  both  plate  and  beam,  the  value  of  a  now  varies  between  bonds. 


6.4  “Boundary”  Conditions 

Because  peridynamic  models  result  in  long  range  forces,  it  is  not  sufficient  to  apply  boundary 
conditions  to  nodes  on  the  relevant  boundary;  nodes  near  the  boundary  must  be  considered  as  well. 


V  Boundary  Region  J 

Figure  6.10:  The  boundary  of  a  peridynamic  model  is  a  region  of  nonzero  thickness 

For  the  peridynamic  beam  we  consider  simple  supports  or  rollers,  and  fixed  or  clamped  sup¬ 
ports.  Simply-supported  beams  are  easy  to  model  because  only  displacement  is  constrained.  To 
add  a  roller  support  to  the  peridynamic  beam,  it  suffices  to  constrain  the  movement  of  the  nearest 
peridynamic  node  in  the  appropriate  degree(s)  of  freedom.  Simulating  a  “clamped”  end  condition 
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is  a  little  less  intuitive.  The  most  basic  way  to  simulate  a  clamped  end  is  to  extend  the  beam  25,  or 
twice  the  horizon,  into  the  clamp.  The  displacement  of  all  of  those  nodes  is  set  to  zero,  or  whatever 
value  is  appropriate  for  a  displaced  or  rotated  clamp.  In  classical  mechanics,  a  clamped  end  can 
be  described  with  a  symmetry  condition,  but  the  two  are  not  peridynamically  equivalent.  Because 
the  classical  beam  is  a  local  model,  material  at  a  clamp  cannot  “see”  distant  material,  so  there  is 
no  way  to  distinguish  between  a  beam  end  that  is  clamped  and  one  that  is  bent  symmetrically  over 
an  appropriate  sawhorse. 

The  loads  we  apply  to  the  peridynamic  beam  include  applied  moments,  point  loads  and  dis¬ 
tributed  loads.  Distributed  forces  may  be  applied  as  expected  to  nodes  in  the  loaded  region.  Point 
forces  may  often  be  applied  directly  to  the  nearest  node,  or  to  the  nodes  immediately  surrounding 
the  point  of  application.  Point  moments  must  also  be  considered  more  carefully  because  the  peri¬ 
dynamic  models  in  this  work,  like  most  peridynamic  models,  do  not  consider  rotational  degrees 
of  freedom  for  peridynamic  nodes.  Rather,  material  rotation  is  the  result  of  the  relative  transla¬ 
tional  displacement  of  multiple  nodes.  It  is  therefore  impossible  to  apply  a  moment  to  a  single 
peridynamic  point.  Instead,  moments  may  be  applied  as  force  couples  to  the  bonds  attached  to  the 
peridynamic  node  nearest  the  location  of  the  desired  moment.  For  example,  if  we  want  to  apply 
a  moment  M  at  point  x,  whose  n  neighboring  points  x'  are  connected  to  x  in  the  undeformed 
configuration  by  bonds  For  an  evenly  discretized  beam,  we  may  distribute  the  moments  by 

Mi  =  M  , 

E 

i=i 

and  apply  them  to  the  corresponding  bonds  by  adding  forces 


Fi  =  Mi  X 


Y{ii) 

Yte)P 


to  each  point  x*  and  subtracting  them  from  the  force  at  x. 

Support  configurations  are  similar  for  two-dimensional  models.  Each  node  along  a  simply- 
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supported  edge  is  constrained  in  one  or  more  directions.  As  with  the  beam  model,  clamped  edges 
are  implemented  by  extending  the  surface  into  the  clamp.  Line  and  pressure  loads  are  treated 
normally. 

6.5  Numerical  Solution  Method 

This  project  uses  Trilinos,  a  collection  of  open  software  libraries,  or  packages,  from  Sandia  Na¬ 
tional  Labs,  including: 

•  Epetra  and  EpetraExt  -  provide  efficient  parallel  data  structures,  particularly  vectors  and 
sparse  matrices 

•  Isorropia  -  provides  load  balancing,  partitioning,  and  matrix  coloring 

•  NOX  -  a  collection  of  large-scale  nonlinear  system  solver  utilities 

•  Py Trilinos  -  a  python  interface  providing  Python  wrappers  for  many  Trilinos  packages,  and 
offering  compatibility  between  numpy.ndarrays  and  Epetra.MultiVectors  [61] 

The  nature  of  discrete  peridynamic  models  results  in  large  numbers  of  parallelizable  compu¬ 
tations.  Efficient  parallelization  is  achieved  using  Epetra  data  structures  for  distributed  variables. 
Model  force  evaluations  are  coded  in  Python,  making  extensive  use  of  the  optimized  routines  in 
the  NumPy  and  SciPy  packages  operating  on  the  distributed  Epetra  objects.  To  obtain  quasistatic 
solutions,  problems  are  coded  into  NOX  objects  and  solved  using  NOX  nonlinear  solvers.  Pre¬ 
liminary  analysis  was  performed  using  a  Newton  Method  solver  on  an  iMac  with  a  3.1GHz  Intel 
Core  i7  processor  and  16GB  RAM,  using  1-4  cores.  Later  work  also  performed  on  Shamu,  a  High- 
Performance  Computing  cluster  at  UTSA,  and  Stampede,  a  High-Performance  Computing  cluster 
at  UT  Austin.  The  nature  of  the  Trilinos  packages  and  the  structure  of  the  code  allow  for  more 
extensive  parallel  computation  without  major  code  changes. 

For  all  but  the  simplest  loading  conditions,  analytical  solutions  to  boundary  condition  problems 
become  complicated.  As  load  conditions  and  material  behavior  become  more  complex,  there  are 
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no  analytical  solutions.  For  comparison,  equivalent  models  are  created  and  analyzed  in  Abaqus 
6.12  to  verify  simple  cases. 
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6.6  Results 


6.6.1  Straight  Beam  Results 

The  simplest  test  case  for  this  model  is  a  linear-elastic  beam  with  a  square  profile.  These  models 
represent  beams  that  are  1cm  thick  and  1cm  wide,  with  a  bulk  modulus  k  =  37.5GPa  and  Poisson’s 
ratio  V  =  Vs.  Each  is  loaded  transversely  with  a  load  of  0.0833N,  except  for  fig.  6.13,  which  shows 
a  beam  loaded  with  a  moment  of  0.0833N-m.  In  beams  simulating  inelastic  behavior,  the  critical 
strain  of  the  material  is  set  to  Ec  =  0.001.  The  elastic  cases  are  compared  to  analytical  solutions 
of  the  Euler  beam  equation  with  appropriate  boundary  conditions.  Even  a  coarse  discretization 
successfully  reproduces  the  shape  of  the  elastically  deformed,  simply- supported  beam  shown  in 
fig.  6.11  deformed  under  uniform  load.  Other  load  types  are  also  possible;  figs.  6.12  and  6.13 


Figure  6.11:  The  uniform-load  elastic  beam  is  accurately  modeled  with  few  nodes 


demonstrate  simply-supported  elastic  beams  with  a  point  load  and  a  point  moment,  respectively. 

It  is  more  difficult  to  accurately  reproduce  the  behavior  of  a  clamped-end  beam.  It  is  evident 
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Figure  6.12:  Simply  Supported  Beam  with  Point  Load 


Figure  6.13:  Simply  Supported  Beam  with  Point  Moment  at  Center 
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from  figs.  6.14  and  6.15  that  the  clamped  end  constraint  requires  far  more  nodes  and  a  smaller 
horizon  to  reproduce  the  results  of  a  classical  elastic  model.  Figure  6.16  shows  a  cantilever  beam 
with  one  clamped  end  and  one  free  end,  deflecting  under  a  uniformly  distributed  load. 


Figure  6.14:  The  clamped  condition  requires  finer  discretization 


As  an  elastic-perfectly-plastic  beam  exceeds  the  elastic  limit  of  its  material,  plastic  zones  begin 


to  grow  on  the  top  and  bottom  of  the  beam’s  cross  section.  This  behavior  is  mimicked  by  the  plas¬ 


ticity  of  the  longest  bond-pairs  described  in  eq.  (5.9),  producing  the  results  shown  in  fig.  6.17.  To 


evaluate  plastic  beam  behavior,  plastic  beam  models  with  identical  material  properties  are  created 


and  analyzed  in  Abaqus  6.12  to  verify  simple  cases.  To  accurately  capture  this  phenomenon  and 


model  beam  plasticity,  a  finer  discretization  is  required. 


A  material  that  is  plastically  deformed  does  not  return  to  its  original  state  when  unloaded.  For 


a  beam  in  bending,  the  residual  deformations  can  be  seen  in  a  beam  that  has  been  loaded  beyond 


the  onset  of  plastic  deformation  and  then  unloaded.  The  Abaqus  model  retains  slightly  more  than 


Vio  of  its  loaded  displacement  after  being  completely  unloaded.  This  result  is  observed  in  the 
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Deflection  under  Uniform  Load 


Figure  6.15:  The  clamped  condition  requires  a  smaller  horizon 


Figure  6.16:  Uniformly  Loaded  Cantilever  Beam 
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Figure  6.17:  The  elastic  perfectly-plastic  beam  requires  finer  discretization 

bond-pair  plasticity  model,  shown  in  figs.  6.18  and  6.19.  Accurate  residual  deformation  modeling 
requires  both  a  relatively  small  horizon  and  a  fairly  large  number  of  nodes.  While  this  makes  for 
a  computationally  expensive  model  in  this  case,  it  allows  for  implementation  of  more  complex 
plasticity  models  without  additional  expense;  a  plasticity  model  that  includes  softening  is  no  more 
effort  to  implement,  and  would  result  in  damage  fields  that  depend  on  the  peridynamic  horizon,  not 
on  the  density  of  the  discretization.  Because  the  peridynamic  horizon  is  ideally  a  material  property, 
this  has  the  effect  of  regularizing  the  solution.  Implementing  a  strain- softening  plasticity  model  in 
a  finite  element  beam  would,  by  contrast,  result  in  the  localization  of  plasticity  into  one  element 
and  depend  strongly  on  the  choice  of  discretization. 

It  is  more  difficult  to  verify  the  brittle  material  model  described  by  eq.  (5.10)  because  brittle 
failure  is  unstable.  When  a  crack  begins,  moment  is  transferred  to  other  bond  pairs,  and  failure 
progresses  until  every  pair  of  bonds  surrounding  a  node  are  broken,  creating  a  hinge  at  that  node. 
This  is  borne  out  by  the  results  in  fig.  6.20,  in  which  “Nodal  Health”  represents  the  fraction  of 
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Figure  6.18:  The  need  for  fine  discretization  is  even  more  apparent  when  representing  residual 
plastic  deformation 

bond-pairs  about  each  node  that  have  never  exceeded  their  critical  angle  and  therefore  have  not 
failed. 

Unlike  a  local  model,  partial  failure  is  observed  at  nodes  near  the  plastic  hinge,  as  pairs  of 
bonds  that  straddle  the  hinge  are  broken.  This  is  an  important  feature  of  peridynamics;  as  in  the 
plastic  beam,  the  damaged  region  of  the  brittle  beam  depends  on  the  peridynamic  horizon  rather 
than  on  the  density  of  the  discretization.  In  a  finite  element  model,  damage  will  occur  between 
elements  or  within  a  single  element.  Either  way,  mesh  refinement  will  eventually  result  in  damage 
localization  to  an  infinitely  small  damage  region. 
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Distance  along  Beam 


Figure  6.19:  Accurately  modeling  residual  plastic  deformation  also  requires  a  small  horizon 
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Figure  6.20:  A  brittle  beam  with  prescribed  center  displacement 
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Node  Health 


6.6.2  Flat  Plate  Results 


The  simplest  test  case  for  the  2D  model  is  a  linear-elastic  Im  by  Im  plate  that  is  simply- supported 
on  all  4  sides  with  a  uniform  transverse  pressure  load  on  the  entire  surface  between  the  supports. 
The  plates  in  this  section  have  a  shear  modulus  k  =  37.5Gpa,  and  unless  otherwise  noted,  a  Pois¬ 
son’s  ratio  of  V  =  1/3.  The  elastic  plates  are  all  loaded  with  a  total  transverse  force  of  937. 5N.  As 
expected  from  an  energy-equivalent  model,  the  slice  along  the  plate’s  centerline  shown  in  fig.  6.21 
demonstrates  good  agreement  between  the  static  deflection  predicted  by  the  bond-pair  model  and 
that  of  classical  linear  elasticity  as  the  horizon  length  shrinks.  This  convergence  only  continues  to  a 
minimum  horizon,  below  which  the  discretized  equation  of  motion  (eq.  (6.2))  ceases  to  accurately 
approximate  the  continuous  integral  formulation  (eqs.  (4.1)  and  (5.1)).  The  minimum  horizon  size 
depends  on  the  discretization;  it  appears  that  three  times  the  node  spacing  is  sufficient,  but  that  a 
horizon  that  is  only  twice  the  node  spacing  is  insufficient.  The  difference  is  evident  in  fig.  6.22, 
which  also  shows  that  results  are  insensitive  to  fineness  of  discretization  once  the  minimum  hori¬ 
zon  criterion  is  met.  Accurate  results  require  a  denser  discretization  than  is  the  case  for  the  elastic 
beams  from  previous  work.  Figure  6.23  illustrates  the  model  converging  to  the  analytical  solution 
as  the  discretization  is  made  finer  and  the  horizon  shrinks. 

The  test  case  for  the  hybrid  model  is  a  similar  simply-supported  square  plate  with  an  additional 
in-plane  tension  load  along  two  opposing  sides.  An  analytical  solution  for  this  combination  of 
uniform  transverse  pressure  and  in-plane  edge  tension  can  be  found  in  Timoshenko’s  book  [74]. 
As  is  mechanically  intuitive,  increasing  in-plane  tension  results  in  decreasing  transverse  displace¬ 
ment,  while  the  opposite  is  true  for  compressive  edge  loading.  Normalized  to  the  the  maximum 
displacement  of  a  transversely-loaded  plate  with  no  in-plane  edge  loads,  the  results  in  fig.  6.24 
show  that  the  hybrid  model  does  a  good  job  of  simulating  the  impact  of  in-plane  tension  on  max¬ 
imum  transverse  deflection.  The  bond-multiple  plate  model  is  motivated  by  the  desire  to  extend 
the  bending  model  to  an  arbitrary  Poisson’s  ratio,  so  the  obvious  test  for  this  model  is  the  same 
as  for  the  bond-pair  model.  When  compared  to  analytical  predictions,  fig.  6.25  demonstrates  the 
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Figure  6.21:  The  Bond-Pair  Model  Converges  on  Accurate  Plate  Deflection  with  Smaller  Horizons 

bond-multiple  model’s  ability  to  simulate  plates  with  Poisson’s  ratios  that  depart  significantly  from 
the  bond-pair  limitation  of  =  ^/s. 

Unlike  the  brittle  beam,  failure  in  a  brittle  plate  need  not  progress  unstably.  To  demonstrate 
the  behavior  of  this  model,  a  controlled-displacement  double-torsion  fracture  test  was  simulated 
with  the  bond-pair  model  using  a  critical  strain  Sc  =  0.001.  A  good  review  of  the  double-torsion 
test  is  available  in  [63].  Figure  6.26  shows  the  setup  of  a  double  torsion  fracture  test.  This  test 
is  particularly  useful  because  it  results  in  a  bending  crack  whose  growth  is  not  unstable.  The  two 
sides  of  the  cracked  plate  act  as  torsion  springs;  as  the  crack  grows  longer,  the  torsion  springs 
grow  longer  as  well,  and  correspondingly  softer.  The  overall  result  is  that,  even  though  the  plate’s 
resistance  to  bending  decreases  as  the  crack  grows,  growth  is  stable  until  the  crack  nears  the  far 
side  of  the  plate.  The  simple  qualitative  results  are  shown  in  fig.  6.27,  colored  by  the  fraction  of 
failed  bond  pairs  around  each  node.  For  each  successive  displacement  load,  the  stable  progression 
of  the  damaged  region  extends  further  into  the  plate.  The  relationship  between  load  displacement 
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Figure  6.22:  Horizon  Must  Include  Sufficient  Nodes 

and  crack  length  is  consistent  with  a  plane-stress  mode  I  fracture  toughness  around  IVGPav^n.  As 
with  the  brittle  beam,  the  region  within  one  horizon  of  the  crack  shows  partial  damage. 

In  the  double-torsion  plate,  the  crack  is  expected  to  progress  straight  across  the  plate.  A  more 
interesting  failure  pattern  can  be  found  if  the  displacement  is  applied  to  only  one  of  the  two  sides 
of  the  pre-crack.  In  a  “single  torsion”  cracked  plate,  the  crack  path  is  expected  to  curve,  and  this 
result  can  be  seen  in  fig.  6.28. 
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Figure  6.23:  The  Bond-Pair  Model  Converges  on  Accurate  Plate  Deflection  with  Finer  Discretiza¬ 
tion 


Figure  6.24:  The  Combined  Model  Accurately  Captures  the  Influence  of  In-Plane  Tension 
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Figure  6.25:  The  Extended  Model  Matches  for  Arbitrary  Poisson’s  Ratio 


, ,  j  Applied  vertical  displacement  load 
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Figure  6.26:  Simple  Double  Torsion  Setup 
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Figure  6.27:  Crack  Progression  in  Double  Torsion  Brittle  Plate 


80 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Damage 


0.005 


1.0 


Initial  Crack 


0.000 


a 

OJ 

a 

(D 

O 

"Hh 

C/5 


-0.010 


-0.015 


Displacement  =  0.020 


0.0  0.5  1.0  1.5  2. 


1.0 

0.8 

0.6 

0.4 

0.2 

^.0 


0.6 


0.4 


0.2 


0.0 


Figure  6.28:  Crack  Progression  in  Single  Torsion  Brittle  Plate 
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Damage 


6.6.3  Irregular  Discretization  Results 


While  it  is  simple  to  discretize  the  shapes  examined  thus  far  so  that  bonds  pair  up  nicely,  that  is 
tougher  to  do  for  most  shapes  we  are  interested  in  analyzing.  To  combine  the  analytical  simplicity 
of  the  previous  shapes  with  a  demonstration  of  the  ability  to  handle  irregular  shapes,  we  discretize 
simple  shapes  in  an  irregular  fashion.  The  first  is  a  return  to  the  simply- supported  beam  under 
uniform  load.  Figure  6.29  shows  that  the  irregularly-discretized  beam  has  the  same  deflection 
under  uniform  load  as  the  regular  discretization. 


Figure  6.29:  Irregular  Beam 


An  example  of  a  regularly-discretized  plate  is  shown  in  fig.  6.30,  while  an  irregularly  dis¬ 
cretized  plate  is  shown  in  fig.  6.31.  To  use  this  plate  in  a  peridynamic  simulation,  we  reduce  each 
element  to  a  single  node  at  the  element’s  centroid  with  volume  equal  to  the  volume  of  the  element. 


82 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Figure  6.30:  The  regular  plate  discretization  allows  pairing  of  opposite  bonds  between  nodes 


Figure  6.31:  The  irregular  plate  discretization  requires  virtual  points  to  form  bond  pairs 
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Figure  6.32:  Virtual  Points  Allow  for  Irregular  Discretizations 
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6.6.4  Curved  Surface  Results 


To  test  the  virtual  point  method  for  curved  surfaces,  we  start  with  the  simple  case  illustrated  in 
fig.  6.33  of  a  ring  in  bending  with  a  beam  cross  section.  Like  the  straight  beam,  the  curved  beam 
is  0.01m  in  both  width  and  thickness,  with  shear  modulus  k  =  37.5Gpa  and  Poisson’s  ratio  u  = 
1/3.  This  Im-diameter  ring  is  subjected  to  a  point  load  of  0.0833N.  A  real-life  example  of  this 
type  of  problem  would  be  a  proving  ring  used  for  force  measurements.  The  plots  in  figs.  6.34 


and  6.35  show  that  convergence  to  the  classical  elastic  solution  is  more  a  function  of  the  shrinking 
horizon  size  than  of  denser  discretization,  though  the  discretization  must  be  dense  enough  to  allow 
sufficient  nodes  in  each  neighborhood. 

Similarly,  the  use  of  virtual  points  in  the  simulation  of  curved  plates  and  shells  can  be  examined 
by  starting  with  the  ring  model  depicted  in  fig.  6.36.  This  ring  is  Im  in  diameter,  0.01m  thick,  and 
0.25m  wide.  The  ring  is  loaded  in  the  same  manner  as  the  beam-based  ring,  but  with  a  line  load 
with  total  magnitude  937.5N.  In  cross-section,  the  ring  looks  like  fig.  6.33. 

Figure  6.37  demonstrates  that,  as  with  the  curved  beam,  the  curved  shell  requires  a  dense 
discretization  to  model.  This  is  mostly  as  a  consequence  of  requiring  a  small  horizon.  Surfaces 
with  less  curvature  are  less  constrained. 
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Figure  6.34:  A  Beam-based  Proving  Ring  in  Varying  Discretizations 
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Figure  6.35:  Horizon  Convergence  in  the  Beam-based  Proving  Ring 
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Figure  6.36:  Discretization  of  a  3D  Proving  Ring 


-  Analytical 

▼  ▼  300  nodes/length,  S  =  0.010 
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Figure  6.37:  Plate-based  Proving  Ring 
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Chapter?:  CONCLUSION 


This  work  develops  what  we  believe  to  be  the  first  state-based  peridynamic  thin  feature  models. 
These  models  allow  for  the  simulation  of  bending  in  peridynamic  beams,  plates,  and  shells.  They 
are  verified  by  comparing  their  strain  energies  to  the  strain  energy  of  classical  models  for  small 
homogenous  deformations,  which  are  matched  for  both  Euler-Bernoulli  beams  and  for  Kirchhoff- 
Love  plates  as  the  peridynamic  horizon  approaches  zero.  In  addition  to  matching  linear  elas¬ 
tic  models  for  these  structures,  it  has  been  shown  that  the  peridynamic  models  also  match  with 
Eringen-style  gradient  elasticity  as  the  horizon  decreases. 

The  bond-pair  version  of  the  bending  model  resists  only  bending  deformation  and  is  limited 
to  plates  with  Poisson’s  ratios  v  =  1/3.  To  simulate  loading  incorporating  forces  both  in-plane 
and  transverse  to  the  beam  or  plate,  the  pure  bending  model  was  combined  with  a  conventional 
peridynamic  model  for  in-plane  deformation.  The  result  is  a  hybrid  bending-extension  model. 
The  Poisson’s  ratio  limitation  was  overcome  by  decomposing  bending  deformation  into  isotropic 
and  deviatoric  components,  allowing  for  a  bond-multiple  isotropic  bending  correction  to  the  de¬ 
formation  energy.  With  this  correction,  the  deformation  energy  was  shown  to  match  that  of  a 
Kirchhoff-Love  plate  for  any  (valid)  Poisson’s  ratio. 

To  further  extend  the  utility  of  the  model,  a  reformulation  was  developed  to  gracefully  include 
curved  and  irregularly  discretized  areas.  By  adding  virtual  points  at  locations  where  no  nodes  is 
present,  the  bond-pair  model  can  be  applied  to  regions  where  bonds  are  not  initially  organized 
into  equal  and  opposite  pairs.  This  includes  curved  surfaces,  where  a  virtual  point  may  be  added 
slightly  above  or  below  the  surface,  as  well  as  irregularly  discretized  regions  in  which  virtual  points 
may  be  added  between  peridynamic  nodes. 

Code  was  written  to  evaluate  both  beam  and  plate  models  for  a  variety  of  support  configu¬ 
rations,  load  configurations,  and  materials.  Although  many  problems  of  interest  result  in  purely 
transverse  deformation,  both  mathematical  models  and  computer  code  are  fully  3D.  The  nature 
of  peridynamic  models  results  in  large  but  easily  parallelizable  problems.  The  code  was  writ- 
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ten  with  this  in  mind,  and  scales  easily  to  large  problems  on  multiple  machines.  Simulations 
run  with  the  developed  models  provide  results  in  agreement  with  conventional  methods  for  sev¬ 
eral  support  configurations,  including  simply-supported,  clamped,  and  free  end  conditions.  Good 
agreement  was  also  reached  for  applied  displacements,  moments,  and  both  point  and  distributed 
loading  conditions.  The  proposed  damage  model  successfully  reproduces  the  impact  of  nonlinear 
elasticity  on  deformation  of  a  rectangular  beam,  and  the  framework  is  laid  to  allow  application 
of  the  same  model  to  a  variety  of  beam  cross  sections.  Brittle  fracture  results  for  both  beam  and 
plate  simulations  are  also  consistent  with  expectation.  A  single-torsion  plate  demonstrates  the 
peridynamic  advantages  of  natural  crack  path  development  without  the  need  to  predict  direction  of 
propagation.  The  hybrid  bending-extension  model  resists  transverse,  in-plane,  and  combined  de¬ 
formation,  allowing  it  to  successfully  reproduce  the  stiffening  effect  of  adding  in-plane  tension  to 
a  transversely-loaded  plate.  With  the  bond-multiple  correction,  the  extended  model  demonstrated 
the  effect  of  varying  Poisson’s  ratio  on  plate  bending.  Because  it  was  inspired  by  a  similar  de¬ 
composition  used  in  extension  models,  there  same  procedure  can  be  used  for  combined  in-plane 
and  transverse  deformations.  Implementation  of  the  virtual  point  extension  greatly  increases  the 
practical  utility  of  the  model,  a  fact  demonstrated  by  accurate  simulations  of  the  displacements  of 
irregularly  discretized  plates  and  beams.  Finally,  the  addition  of  virtual  points  allows  application 
to  curved  features,  as  demonstrated  by  successful  simulation  of  beam  and  plate  based  rings. 

The  state-based  nature  of  these  models  allows  them  to  resist  bending  as  a  basic  deformation 
mode,  in  a  way  that  is  fundamentally  different  from  previous  peridynamic  models,  both  3D  solid 
and  thin-feature.  While  the  examples  in  this  work  focus  on  conventional  elastic  and  inelastic 
behavior,  treating  bending  as  a  basic  deformation  opens  the  possibility  of  modeling  thin  features 
whose  resistance  to  bending  is  not  simply  a  function  of  stresses  that  vary  through  the  thickness 
to  produce  a  moment,  and  suggests  a  wide  variety  of  possible  applications:  Graphene  sheets  are 
only  one  atom  thick,  so  their  bending  behavior  cannot  be  the  result  of  through-thickness  variation 
in  stress.  Some  biological  membranes  resist  isotropic  bending  but  not  deviatoric  stresses.  Highly- 
structured  metamaterials  may  also  have  bending  properties  that  cannot  be  derived  solely  from  bulk 
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solid  properties. 


Peridynamic  modeling  is  a  technique  that  is  swiftly  growing  in  popularity,  especially  among 
those  concerned  with  the  propagation  of  material  failure.  The  bending  behavior  of  thin  features  is 
a  critical  part  of  many  engineering  analyses  that  has  been  largely  out  of  reach  for  reasonably- sized 
peridynamic  models.  These  peridynamic  bending  models  extend  the  domain  of  problems  to  which 
peridynamic  modeling  can  be  applied. 

There  are  many  directions  in  which  the  models  developed  in  this  work  can  be  improved  upon. 
While  the  results  presented  here  are  both  illustrative  and  persuasive,  they  do  not  represent  a  mathe¬ 
matically  rigorous  demonstration  of  the  convergence  properties  of  the  bending  models,  and  future 
work  focusing  specifically  on  discretization  and  convergence  issues  would  be  valuable.  Both  the 
plastic  and  brittle  material  models  developed  in  this  work  are  proofs  of  concepts,  demonstrating  the 
capabilities  of  peridynamics  and  the  new  bending  models.  It  will  be  necessary  to  develop  a  thermo¬ 
dynamically  consistent  failure  model  before  using  this  model  to  make  predictions  regarding  failure 
propagation,  especially  when  rapid  progression  is  expected.  Ideally,  the  material  parameters  of 
such  a  model  could  be  determined  from  existing  datasets  and  material  parameters. 

Because  directional  composites  are  a  popular  material  choice  for  thin  features,  it  would  be 
worthwhile  to  extend  the  bending  models  to  simulate  non-isotropic  plates.  Other  useful  avenues 
of  inquiry  include  combining  these  2D  models  with  peridynamic  solid  models  to  analyze  features 
comprised  of  both  thick  and  thin  parts. 
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Appendix  A:  FRECHET  DERIVATIVE 


A.l  Definition 

The  derivative  of  a  function  of  a  state  is  defined  by  Silling  in  [66]  as  follows: 

Let  be  a  function  of  a  state,  ^'(•)  :  Am  ^n-  Suppose  there  exists  a  state-valued 
function  denoted  G  Am+n  such  that  for  any  A  G  Am  and  any  A  A  G  Am, 

^(A  + AA)  =  ^(A)  +  V^(A)  •  AA  +  o(||AA||).  (A.l) 

Then  ^  is  said  to  be  differentiable  and  is  called  the  Frechet  derivative  of 

This  is  a  fairly  straightforward  way  of  defining  a  derivative  with  respect  to  a  state.  Because  the 
force  vector-state  and  deformation  vector-state  are  work  conjugate,  the  force  vector-state  can  be 
determined  by  taking  the  Frechet  derivative  of  energy  with  respect  to  the  deformation  vector-state. 

A.2  Bond-Pair  Force 

For  the  bond-pair  model,  we  derive  the  bond  force  function  from  the  bond-pair  energy  function 

T®  =V«,(Y(e)) 

w  =  u.(«) a  [l  +  cos(«(y{«>,  Y{-0))' 

»(y{«)  +  AY{«))  =u>(«)  tt  1^1  -f  COS  (o(y{«>  +  ay{«>,y{-o))] 

Vro(Y(«))  .  AY(|)  =  u>(y{?)  +  AY{?))  -  u>(y{?)) 

=  U.K)  a sm(o(Y(«>.  Y{-«>) )  [o(y(«>  +  AY(«>.  Y{-«>)  -  o(y{0,  Y{-0)' 

[e(Y«>  +  AY«>,Y{-«))  -0(y(«),Y{-O)]  =  ^^•«(y{«),Y(-«>) 
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To  determine  the  6  direction  vector,  we  must  construct  a  vector  that  is  normal  to  Y(^)  and  that 
is  in  the  plane  containing  both  Y(^)  and  Y(— ^).  The  cross  product  of  Y(^)  and  Y(— is 
a  vector  normal  to  that  plane,  so  any  vector  normal  to  that  cross  product  will  be  in  the  correct 
plane.  Therefore,  the  vector  Y(^)  x  [Y(^)  x  Y(— ^)]  is  both  normal  to  Y(^)  and  is  in  the  plane 
containing  both  Y(^)  and  Y(— ^).  Normalizing  gives  us  the  9  direction  vector: 


»(y{«>.Y(-«>) 


Y{«>  X  [y(«)  X  Y(-(,) 

Y{i)  1 1  Y«)  1 1  Y(-?)  I  sill  («  (Y(e),  Y{-« ))  ) 


We  combine  all  of  these  to  get  the  expression  for  bond  force  found  in  eq.  (5.1). 


T(^) 


^(0 


Y(^) 

|Y(^)I|Y(^)I 


Y(0  Y(-0 

Y(^)|  |Y(-0I 


A.3  Isotropic  Bending  Correction 


To  derive  the  bending  “pressure”  force,  we  start  with  the  isotropic  energy  discrepancy 


TTX*  n  3l^  —  1  _2 


with 


k(Y)  =  - 


m 


m 


0  Jo 


e 


S  /‘27r  'Y(^) 


0  Jo 


e 
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Because  k  is  itself  a  vector-state,  we  will  need  to  begin  with  the  change  in  n  with  respect  to  Y  and 
carry  the  result  through  to  find  the  change  in  W*. 


R,{Y  +  AY) 


2 
m 


«(Y)  +  - 
m 


0  JO 


e 


lY*  ( Y  +  AY) 


V1Y*(Y) 


2/i 


12  1  -  z/ 
Su  —  1 
12  1-u 


ft(Y)  ^  «(Y) 


(*6  p27t 


+- 


lY*  (Y)  +  2/1 


JO  JO 

2>v  —  1  4  u{^)  _ 
12  1  —  i'  m 


^(Y)  .  AY  +  o(||AY||) 


8fih^3u  -luj{^) 

=  m 


This  demonstrates  the  bond-length  dependent  “pressure”  applied  to  each  point  in  the  neighborhood 
of  a  point  with  average  curvature  k. 
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Appendix  B:  NOTATION 


Peridynamics  is  a  new  field,  and  different  authors  use  a  variety  of  notations  to  represent  a  variety 
of  concepts.  Although  chosen  to  be  as  consistent  as  possible  with  other  authors,  the  following 
notation  is  therefore  not  portable  to  other  works. 

Table  B.l:  Peridynamic  Notation 


Notation 

A 

b 

b 

c 

D 

E 

e 

E 

F 

F 

F 

f 

r 

n 

I 

K 

k 

K 


Meaning 
area 
beam  width 
body  force  vector 
bond  stiffness 
plate  flexural  rigidity 
Young’s  modulus 
bond  stretch 
isotropic  portion  of  bond  stretch 
deviatoric  portion  of  bond  stretch 
applied  force 
deformation  gradient  tensor 
approximated  deformation  gradient  tensor 

force  vector 
isotropic  bending  “pressure” 
neighborhood 
second  moment  of  area 
Eringen’s  nonlocal  modulus 
bulk  modulus 
shape  tensor 
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Table  B.l:  Continued 


Notation 

I 

M 

m 

O 

P 

P 

P 

P,q,r 

Q 

s 

t 

t 

tlD 

T 

u 

V 

V 

w 

w 

W 

W* 

x,p,z 

x,q 


Meaning 

beam  length  represented  by  node  (ID) 

moment 

peridynamic  normalization  constant 
order  notation 
pressure 

distributed  pressure  load  (2D) 
first  Piola-Kirchhoff  stress  tensor 
deformed  bond  vectors 
distributed  load  (ID) 
bond  stretch 
thickness 
Eringen  nonlocal  stress 
ID  Eringen  nonlocal  stress 
force  vector  state 
correction  force  vector  state 
displacement  vector 
shear  force 
displacement  in  y  direction 
displacement  in  2  direction 
bond  energy,  bond-pair  energy 
strain  or  deformation  energy  density 
energy  discrepancy 
coordinate  vectors 
undeformed  location  vectors 
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Table  B.l:  Continued 


Notation 


Meaning 


y 

Y 

OL 

5 

Ax 

e 

€ 


e 

Oc 

(^3D,  O2D 


K 

R 

h 

d'^w 

l^xx 

=  Kl  = 

aw 

hvy  = 

Kyy 

=  K,2  = 

dy"^ 

d'^w 

l^xy  ~ 

=  ^3 

dxdy 

V 


ix 


ii 


deformed  location  vector 
deformation  vector  state 
peridynamic  bending  stiffness 
isotropic  bending  stiffness 
horizon 

distance  between  nodes  in  regular  discretization 

uniaxial  strain 
strain  tensor 
critical  strain 
deviatoric  strain 
bond-pair  angle 
critical  angle 
dilation 
beam  curvature 
nonlocal  curvature 
vector  isotropic  curvature 

Plate  Curvatures 


shear  modulus 
Poisson’s  ratio 
undeformed  bond  vectors 
i-th  of  n  vectors  to 
z-th  component  of  vector  ^ 
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Table  B.l:  Continued 


Notation  Meaning 


e 

magnitude  of  vector  ^ 

p 

density 

p 

radius  of  curvature 

O'lD 

ID  local  stress 

a 

Cauchy  stress  tensor 

(r/) 

scale  parameter 

0 

beam  deflection  angle 

0 

bond  orientation  angle 

Q 

undeformed  body 

Q 

classical  strain  energy  density  function 

UJ 

weight  function 
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